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We have developed an efficient simulation tool ’GOLLUM’ for the computation of electrical, spin and ther¬ 
mal transport characteristics of complex nanostructures. The new multi-scale, multi-terminal tool addresses a 
number of new challenges and functionalities that have emerged in nanoscale-scale transport over the past few 
years. To illustrate the flexibility and functionality of GOLLUM, we present a range of demonstrator calcula¬ 
tions encompassing charge, spin and thermal transport, corrections to density functional theory such as LDA+U 
and spectral adjustments, transport in the presence of non-collinear magnetism, the quantum-Hall effect, Kondo 
and Coulomb blockade effects, finite-voltage transport, multi-terminal transport, quantum pumps, supercon¬ 
ducting nanostructures, environmental effects and pulling curves and conductance histograms for mechanically - 
controlled-break-j unction experiments. 

PACS numbers: 


I. INTRODUCTION 

The development of multi-functional codes capable of pre¬ 
dicting quantum transport properties of complex systems is an 
increasingly active field of research 1 5 . This is driven in part 
by the top-down scaling of the active elements within CMOS 
devices, for which quantum effects are becoming important. 
It is also driven by the bottom-up demands of communities 
working on single-molecule electronics and low-dimensional 
systems, where structures and molecules of increasing size 
and complexity are of interest. In particular, multi-functional 
codes are needed to describe the fundamental properties of 
quasi-two-dimensional materials such as graphene, silicene, 
germanene and their integration into workable devices. The 
need to understand the interplay between all of the above 
structures and their surrounding environments creates further 
demands for such codes. At a more fundamental level, over 
the past forty years, a ’standard model’ of electron transport 
has been developed, based on computing the scattering matrix 
of quantum systems connected to external sources and there 
is a need for a universal code which describes the many real¬ 
izations of such systems under a common umbrella. 

A key task of any quantum transport code is to start from the 
Hamiltonian describing a system and calculate the quantum- 
mechanical scattering matrix S, from which a wide range of 
measurable quantities can be predicted. Unfortunately for 
most nanoscale systems of interest, the full many-body atom¬ 
istic Hamiltonian is too complex to allow this task to be com¬ 
pleted and therefore one usually resorts to a description based 
on a mean-field Hamiltonian. The system of interest is then 
composed of a scattering region, connected to external crys¬ 
talline leads, which are in turn connected to external reser¬ 
voirs. The problem of computing the scattering matrix for 
such a system described by an arbitrary mean-field Hamil¬ 
tonian is solved in reference. The main question therefore 
is how to obtain the correct mean-field Hamiltonian. The 


simplest mean-field approach to describing quantum trans¬ 
port through nanostructures is to build a tight-binding Hamil¬ 
tonian, which reproduces key electronic properties near the 
Fermi energy. This approach has been available for more 
than half a century and is still popular today when describ¬ 
ing generic properties of materials such as graphene 7 . 

Tight binding parameters can be obtained by fitting to 
known band structures and then varied spatially to describe 
external fields and other perturbations. However such an ap¬ 
proach does not easily capture the effects of interfaces be¬ 
tween different materials or edge terminations of finite-size 
systems, whose properties are distinct from those of bulk ma¬ 
terials. Nor does it easily describe finite-voltage effects. To 
capture these additional features of inhomogeneous systems, 
a more material-specific approach is needed. This problem 
was solved in part by the non-equilibrium Green’s function 
technique 8 24 , that combines with density functional theory 
(DFTjMEO obtain the self-consistent mean-field Hamilto¬ 
nian of the system subject to a finite bias voltage and from 
it, the lesser Green’s functions providing the non-equilibrium 
electronic density and current. This approach is utilized 
within the SMEAGOL cod^ 9 , which was the first to describe 
spin-dependent and finite-voltage transport properties of sys¬ 
tems with inhomogeneous magnetic moments and in the pres¬ 
ence of spin-orbit scattering. 

It is almost 10 years since the release of SMEAGOL 
and during this period, we have developed a new code 
with increased speed, versatility and functionality, which is 
particularly suited to the modeling of larger-scale nanos¬ 
tructures, interacting with their environments. This new 
code is called GOLLUM and will be freely available 
from http://www.physics.lancs.ac.uk/gollum within the com¬ 
ing weeks. Our previous experience in the development of 
SMEAGOL 9 10 allowed us to understand that non-equilibrium 
transport codes are quite difficult to handle, in part because of 
their complex input data structures, which can create a steep 
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learning curve, and also because they carry very heavy com¬ 
putational demands. As a consequence, we have devised the 
new code GOLLUM to be more user friendly, with simple and 
easy to understand input and output structures, and having no 
accuracy parameters to tune. We present now a short sum¬ 
mary of the features and functionalities of the two programs 
to better appreciate its differences. 

SMEAGOL is a NEGF program that computes the charge 
and spin transport properties of two-terminal junctions sub¬ 
ject to a finite voltage bias. SMEAGOL cannot read a user- 
defined tight-binding Hamiltonian. Instead, it reads the mean- 
field Hamiltonian from the program SIESTA 27 and is tightly 
bound to the old versions of it. SMEAGOL can read from 
SIESTA Hamiltonians carrying non-collinear spin arrange¬ 
ments as well as the spin-orbit interaction. SIESTA and 
SMEAGOL have indeed been used successfully to simulate 
the magnetic anisotropies of atomic clusters 28 30 and the spin 
transport functionalities of several atomic chains and molec¬ 
ular junctions subjected to strong spin-orbit interaction 3132 . 
However, SMEAGOL does not profit from other recent den¬ 
sity functionals. Examples are the van der Waals family of 
functionals or those based on the LDA+U approach. 

GOLLUM is a program that computes the charge and spin, 
and the electronic contribution to the thermal transport prop¬ 
erties of multi-terminal junctions. In contrast to NEGF codes, 
GOLLUM is based on equilibrium transport theory, which 
means that it has a simpler structure, it is faster and con¬ 
sumes less memory. The program has been designed for user- 
friendliness and takes a considerable leap towards the realiza¬ 
tion of ab initio multi-scale simulations of conventional and 
more sophisticated transport functionalities. 

The simpler interface of GOLLUM allows it to read model 
tight-binding Hamiltonians. Furthermore, GOLLUM has 
been designed to interface easily with any DFT code that 
uses a localized basis set. It currently reads information 
from all the latest public flavors of the codes SIESTA 27 and 
FIREBALL®. These include functionals that handle the spin- 
orbit or the van der Waals interactions, or that include strong 
correlations in the spirit of the LDA+U approach. Plans 
to generate interfaces to other codes are underway. Two- 
and three-dimensional topological materials display fascinat¬ 
ing spin transport properties. GOLLUM can simulate junc¬ 
tions made of these materials either using parametrized tight- 
binding Hamiltonians 34 35 , or DFT®. 

DFT does not handle correctly strong electronic correla¬ 
tion effects, that are inherent many nano-scale electrical junc¬ 
tions. As a consequence, a number of NEGF programs like 
SMEAGOL underestimate such effects. GOLLUM includes 
several tools to handle strong correlations. These include the 
above-mentioned interface to the versions of SIESTA contain¬ 
ing the LDA+U functional. A second tool uses a phenomeno¬ 
logical but effective approach called the scissors correction 
scheme. A third tool maps the DFT Hamiltonian into an 
Anderson-like Hamiltonian that is handled with an impurity 
solver in the spirit of dynamical mean field theory. 

The lighter computational demands required by GOLLUM 
make it possible to construct conductance statistics relevant 
to break-junction and STM measurements of single-molecule 


conductances, therefore making closer contact with experi¬ 
ments. GOLLUM also incorporates an interface with some 
classical molecular dynamics programs, which enables it to 
handle interactions with the environment. 

GOLLUM makes use of the concept of virtual leads , that 
allows it to integrate easily a wide range of phenomena by the 
use of tight-binding Hamiltonians. These include spintronics, 
superconductivity, Kondo physics and topological phases. In 
contrast with SMEAGOL, which only computes the magni¬ 
tudes of transmission coefficients of two-terminal junctions, 
GOLLUM has access to the full scattering matrix of a multi¬ 
terminal junction, enabling it to compute scattering ampli¬ 
tudes, phases and Wigner delay times and thereby describe 
the properties of quantum pumps. 

Even though GOLLUM is based on equilibrium trans¬ 
port theory, our experience with the use of the NEGF code 
SMEAGOL has enabled us to incorporate non-equilibrium 
physics into the mean-field Hamiltonian. GOLLUM has 
therefore the ability to compute non-equilibrium current- 
voltage curves. 

In this article, our aim is to describe the structure of the 
code and then present a set of demonstrator calculations. The 
latter will illustrate the additional functionality and versatil¬ 
ity of GOLLUM and at the same time constitute a set of new 
results for the transport properties of selected structures. All 
of the functionalities that will be discussed below are avail¬ 
able either in the current public version of the code, or in the 
current development version, that will be made public in the 
autumn of 2014. 

The layout of this article is as follows. In Section II, we de¬ 
scribe the theoretical approach behind the program and outline 
the theoretical and practical details of the current implemen¬ 
tation. The section starts with a detailed description of the 
generic two-probe and multi-probe junction setups available 
within GOLLUM and introduces the terminology that will be 
used throughout the article. This is followed by subsections 
describing the determination of the surface Green’s function 
of each current-carrying lead and the full scattering matrix. 
We then introduce a convenient method that allows us to de¬ 
scribe finite-voltage non-equilibrium effects. A subsection ex¬ 
plaining the concept of virtual leads enables us to describe hy¬ 
brid structures containing non-collinear magnetism or super¬ 
conductivity within the scattering region. Two additional sub¬ 
sections explain two facilities included in GOLLUM that en¬ 
able us to describe electronic correlation effects beyond DFT, 
including Coulomb blockade and Kondo physics. We then 
show how to include a gauge field in the GOLLUM Hamilto¬ 
nian. A final subsection explains the multi-scale methodology 
used to describe large-scale junctions and environmental ef¬ 
fects, using a combination of classical molecular dynamics for 
the environment and quantum transport for the central scatter¬ 
ing region. 

In section III, we present the details and results of the 
simulations of a series of sixteen different demonstrator sys¬ 
tems. The purpose of each demonstrator is to present one 
or more of the functionalities of GOLLUM. We start with 
a few model junctions, described by tight-binding Hamilto¬ 
nians, which show basic capabilities and demonstrate how 
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easily and flexibly the program can analyze non-trivial phys¬ 
ical effects. These include two and four-terminal normal- 
metal junctions, a two-dimensional system showing the quan¬ 
tum Hall effect, and two hybrid structures, the first contain¬ 
ing two superconducting islands sandwiched by two normal 
metal electrodes, and the second containing a non-collinear 
spin structure. There follow a series of DFT-based calcula¬ 
tions, that describe spin-active junctions; graphene junctions, 
where the use of van der Waals functionals is crucial and junc¬ 
tions enclosing metallo-organic molecules, where a treatment 
of strong-correlations beyond DFT is mandatory. These are 
handled using three different approaches. In the first, the prop¬ 
erties of metalloporphyrin junctions are described using the 
LDA+U methodology; in the second, the LDA spectrum of 
an OPE derivative is adjusted to improve the agreement with 
experimental data; in the third, we describe Coulomb block¬ 
ade and Kondo features of model or simple gold junctions. 
The next demonstrator shows how GOLLUM can compute the 
thermoelectric transport properties of a junction.This is fol¬ 
lowed by a discussion of the transport properties of a carbon- 
nanotube-based four-probe junction. The next three exam¬ 
ples require the use of multi-scale techniques, where we use 
a three-step methodology described later in the article. The 
first demonstrator describes how liquid environmental effects 
modify the transport properties of a single-molecule junc¬ 
tion. The second example demonstrates that single strands of 
DNA can be trans-located though graphene nanopores, where 
the strands effectively gate the nanopore structure yielding a 
highly sensitive DNA sensor based on field-effect-transistor 
concepts. A final demonstrator shows how GOLLUM can 
compute full sequences of pulling and pushing cycles in 
single-molecule junctions resembling the opening and closing 
cycles of Mechanically Controllable Break Junction (MCBJ) 
experiments, enabling the construction of theoretical conduc¬ 
tance histograms. The final demonstrator shows how GOL¬ 
LUM has access to the phase of the full scattering matrix, and 
describes non-trivial quantum pumping effects related to the 
phase evolution of the scattered wave function. A concluding 
section summarizes the features delivered by the program and 
an appendix illustrates our method to decimate Hamiltonians 
and overlap matrices. 


II. THEORETICAL APPROACH 
A. Description of the transport methodology 

1. The generic setup and construction of the Hamiltonian 

GOLLUM describes open systems comprising an extended 
scattering region (colored dark blue in Figs. 1 and 2) con¬ 
nected to external crystalline leads (colored light blue in Figs. 
1 and 2). Depending on the problem of interest and the lan¬ 
guage used to describe the system, the material (M) of interest 
forming the central part of the scattering region could com¬ 
prise a single molecule, a quantum dot, a mesoscopic cavity, a 
carbon nanotube, a two-dimensional mono- or multi-layered 
material, a magneto-resistive element or a region containing 
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FIG. 1: (Color online) (a) Schematic plot of a four-terminal device, 
which includes an extended scattering region and four leads. 


one or more superconductors. 

Figure 1 shows an example of a 4-lead system whose cen¬ 
tral scattering region (generically labelled M throughout the 
paper) is a molecule. It is important to note that in an ac¬ 
curate ab initio description of such a structure, the proper¬ 
ties of the leads closest to the molecule (or more generally 
the central scattering material) will be modified by the pres¬ 
ence of the central scattering region (M) and by the fact that 
the leads terminate. In what follows, we refer to those af¬ 
fected portions of the leads closest to the central scatterer as 
’branches’ and include them as part of the ’extended scatterer’ 
(denoted EM throughout this article). Consequently within 
GOLLUM, a typical structure consists of an extended scat¬ 
terer (EM), formed from both the central scatterer (M) and the 
branches. The extended scattering region is connected to crys¬ 
talline current-carrying leads of constant cross-section, shown 
in light blue in the Figs. 1 and 2. For an accurate description 
of a given system, the branches are chosen to be long enough 
such that they join smoothly with the (light blue) crystalline 
leads. Crucially, the properties of this interface region be¬ 
tween the central scatterer M and the leads are determined by 
their mutual interaction and are not properties of either M or 
the electrodes alone. 

Fig. 2 shows a two-terminal device in more detail and in¬ 
troduces further terminology to be used throughout the paper. 
The regions in light blue are called electrodes or leads and are 
described by perfect periodic Hamiltonians subject to chosen 
chemical potentials. Each lead i is formed by a semi-infinite 
series of identical layers of constant cross section, which we 
refer to as principal layers (PLs). Fig. 2 shows only two PLs 
per lead (colored white), although an infinite number is im¬ 
plied. Furthermore in the figure, the leads are identical and 
therefore the lead index i has been dropped. These PLs are de¬ 
scribed mathematically by intra-layer Hamiltonians Hq. PLs 
must be chosen so that they are coupled only to their near¬ 
est neighbors by the Hamiltonians H{, which means that in 
the presence of long-range couplings, a PL may contain more 
than one longitudinal unit cell of the lead. Then, if each PL 
contains N l orbitals, then Hq and H{ are square N z x N l ma¬ 
trices. The extended scatterer (EM) in dark blue is composed 
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FIG. 2: (Color online) (Top) Schematic two-terminal device, where 
electrons are driven from the left to the right lead through the Ex¬ 
tended scattering region. The leads are possibly kept at chemical po¬ 
tentials / il,r = =be V/2, where V is an applied bias. (Bottom) Each 
lead is composed of an infinite chain of identical PL with Hamilto¬ 
nian Ho coupled with each other via coupling Hamiltonians Hi . The 
extended scattering region comprises the actual central scattering re¬ 
gion and several PL in each branch up to the TPL. The TPL connect 
the EM region to the leads. The central scattering region consists in 
this example of the electrodes surfaces and a molecule. 


of a central scattering region (M) and branches. Each branch 
contain several PLs. These PLs have an identical atomic ar¬ 
rangement as the PLs in the leads. However, their Hamilto¬ 
nians differ from Hq and H\ due to the presence of the cen¬ 
tral scattering region. PL numbering at each branch starts at 
the PL beside the central scattering region. The outermost 
PL at each branch of the EM region is called the terminating 
principal layer (TPL) and must be described by Hamiltoni¬ 
ans Hq TPL and H l { TPL which are close enough to Hq and 
H\, to match smoothly with the corresponding lead Hamilto¬ 
nian. For this reason, GOLLUM requires that the EM con¬ 
tain at the very least one PL. The central scatterer (M) itself 
is described by an intra-scatterer Hamiltonian H ^ and cou¬ 
pling matrices to the closest PLs of the branches. In the ex¬ 
ample in Fig. 2, the central scattering region M comprises a 
molecule and the atoms forming the electrode surfaces. The 
surfaces in GOLLUM include all atoms belonging to the elec¬ 
trodes whose atomic arrangements cannot be cast exactly as a 
PL, due to surface reconstructions, etc. For simplicity, Fig. 2 
shows the case of a symmetric system, although no such sym¬ 
metries are imposed by GOLLUM. All Hamiltonians are spin- 
dependent, but again for notational simplicity, the spin index 
<7 will not be written explicitly here. 

This means that the Hamiltonian W for a given lead i can 


be written as the semi-infinite matrix: 
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When using a non-orthogonal basis set, overlap matrices must 
be defined with the same structure as the Hamiltonian matri¬ 
ces: S'o ±i an d It is convenient to introduce the notation 

Kl ± i = H^ ±1 -ESl ±1 

V = n i ~ES i (2) 

For notational simplicity, from now we will consider the case 
where all leads are equal so that the super-index i can be omit¬ 
ted, although GOLLUM imposes no such restrictions. The 
program can assume either open or periodic boundary condi¬ 
tions in the plane perpendicular to the transport direction. In 
this last case, unit cells are chosen in the plane perpendicular 
to the transport direction and the Hamiltonians and overlap 
matrices acquire a specific dependence on the transverse k- 
vector, k±. 


K;h(k±) = E K oii(R±) e ik± ' R± (3) 

R± 


where /i and v label the N orbitals in the unit cell (ie PL) 
at the origin of R±, while z/ denote orbitals equivalent to z/, 
but placed in adjacent unit cells located at transverse positions 
R±. For Kq, /i and v must belong to the same PL n, while 
for Ki, z/ must belong to the PL n' = n + 1. Finally R± 
are vectors in the two-dimensional Bravais lattice, joining the 
unit cell taken as origin with its neighboring unit cells. 

To illustrate how an EM is connected to leads, we now con¬ 
sider the 4-lead example of Fig.l, where we assume that the 
TPL is the third PL in each branch. To describe such a multi¬ 
terminal setup, the Hamiltonian matrix /C EM of the EM is 
arranged in a non-conventional way. The first matrix block 
corresponds to the central scattering region K the second 
matrix block corresponds to the PLs in the EM branch con¬ 
necting to lead 1, the third matrix block to the PLs in the EM 
branch connecting to lead 2, and so on. 
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Finally, the EM described by JC EM and the leads described leads, 
by JC l are coupled via matrices JC lM to yield a four-terminal The transport properties of the junction are encapsulated in 
junction described by the infinite matrix: its scattering matrix S, which can be obtained by computing 
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( 6 ) 


This arrangement of the Hamiltonian enables the straightfor¬ 
ward generalization of the approach to an arbitrary number of 


by solving the infinite system of equations 

-jcg = x ( 8 ) 

This equation can be simplified by replacing the semi-infinite 
lead Greens functions Q x by their surface Green’s functions 
G x s 0 , whose dimensions are N x N. The remaining system 
of equations takes the form 
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where the coupling Hamiltonians are now 

k 1m = (o oo o o o|o o iG 1 0 0 0 1 0 0 0 1 0 0 «jo) 

and the surface Green’s functions of the isolated leads G l s 0 
can be obtained as described in Section II.A.3 below. The 
equation for the full Green’s function can be written in a more 
compact form as 


/ G,90 

I _^coup^t —JC EM 


I G s G SM 

I qMS gEM 


/. ( 11 ) 


The scattering matrix can be computed using the Green’s 
functions matrix elements connecting the different leads, 
which can be obtained by inverting only the upper matrix box 
in Eq. ([9]) 

G s = (Ggfi + K coup /C EM (iT X)up ) 1 j _1 (12) 

In contrast, access to the local electronic and current densities 
at the EM region is obtained from 

g EM = - (/C EM + (K coup ) t Gg^ 0 K coup ^j 1 (13) 

The above expressions for the Hamiltonians are very gen¬ 
eral. Any appropriate tight-binding Hamiltonian could be in¬ 
troduced by hand to allow computation of the transport prop¬ 
erties of a parametrized model. Alternatively, any DFT code 
using localized basis sets can provide them. In this case the 
DFT program produces the Hamiltonians and Fermi energy of 
the EM region H EM , S EM and E EM , and of each lead Hq ±1 , 
Sq ±1 and E l F in separate runs. GOEEUM has an interface 
to the latest versions of the DFT program SIESTA (SIESTA 
3.1, SIESTA VDW and SIESTA LDA+U) and of FIREBALL 
and more interfaces will be developed in the future. Spin de¬ 
grees of freedom in spin-active systems are handled as fol¬ 
lows: if the spins are all collinear, then we compute sepa¬ 
rate Hamiltonians and perform separate transport calculations 
for the spin-up and -down degrees of freedom. However, if 
the junction has non-collinear spins or is subject to spin-orbit 
interactions 2 ®!, then the spin-up and down-degrees of free¬ 
dom are regarded as two distinct sets of orbitals in the Hamil¬ 
tonian, whose distinct labels allow the computation of spin- 
currents or magneto-resistive behaviors. 


Unit cell 

Hj Hj Hj Hj Hj Hj Hj Hj 
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FIG. 3: (Color online) Infinite system used to generate the leads 
Hamiltonians Kq and K \. A positive direction is defined to be to¬ 
wards the scatterer. 


2. Generating the lead surface Green functions G 5 0 

Each of the lead Green’s functions G l s 0 is determined fol¬ 
lowing the procedures described in Refs. 16191371 with some 
minor modifications. To do so, we start by associating to each 
semi-infinite lead i a periodic infinite system, whose unit cell 
contains a single PL, as sketched in Fig. ([3]). Kq and K\ 
can then be created for this infinite system by hand as model 
Hamiltonians, or can be generated by a DFT program in a ded¬ 
icated simulation. Notice that we will drop the i super-index 
until the end of the section for simplicity. 

By expanding the Bloch eigenstates of the infinite system 
in a localized basis set 

\*{k)) = Y j e ikn c,{k) \<p(n,ri) (14) 


where n, p are indices for its unit cells and the orbitals within 
them, and k is a dimensionless, longitudinal Bloch wave- 
vector, the following secular N x N equation can be deduced 

(K 0 + K ie ik + K_ ie ~ ik ) C{k) = 0 (15) 

where the column vector C(k) contains the wave-function 
coefficients c p (k). The above equation is usually solved by 
choosing a wave vector k and solving for the eigen-energies 
and the corresponding eigenvectors. However, in the present 
transport problem, we do the opposite: we choose the energy 
E and solve for the allowed wave vectors and corresponding 
eigenstates. For a given energy E , the above equation has 
2N solutions with either real or complex wave vectors k p , 
p = 1,..., 27V. To obtain the latter, the equation can be recast 
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where In and Oat are the N x TV identity and zero matrices, 
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and C(k p ) is a N-component column vector. GOLLUM 
solves equation ( [T6| ) as it is superior in numerical terms to 
equation We compute the group velocities of the states 
corresponding to real wave vectors as 


v(k p ) 


(*(MN*(M) r 

(y(k p )\*(k p )) 

c{k p y { R i eikp - K -i e ~ ikp ) c iK ) 

* C(k p ) t (5o + sl e ifc * + 5_1 e- <fc p) C(k p ) 


Note that has units of energy and therefore the real, 

fully-dimensioned group velocity is v(k p )(a/h ), where a is 
the spacing between neighboring PLs in a given lead. 

We now divide the 2N wave vectors obtained from Eq. 
( [16] ) into two sets, each containing N values.The first set de¬ 
noted {k p } are real (complex) wave vectors that have v p > 0 
(Im(fep) > 0) and therefore propagate propagate (decay) to 
the right of the figure. They are consequently called positive 
open (closed) channels. The second set denoted {k p } are real 
(complex) wave-vectors that have v p < 0 (Im(fc p ) < 0) prop¬ 
agate (decay) to the left of the figure. They are called negative 
open (closed) channels. 

We introduce the dual vectors D(k p ), D(k p ), which satisfy 
D(k P Y • C(k q ) = S Piq and D(k P Y • C(k q ) = 5 p , q . These can 
be found by inverting the N x N matrices 


Q = (C(ki), ...,C(k N )) = (Ci,..., Cat) 

Q = (C(k 1 ),...,C(k N )) = (C 1 ,...,C N ) 
(D 1 ,...,D n ) = (Q- 1 ) t 

(D u ...,D n ) = (Q- 1 ) t (19) 

The above eigenvectors can be used to construct the following 
transfer matrices 


N 

T = J 2 C n e ik -Di 

1 

N 

T = Yj C n e~ ikn D\ 

1 

These transfer matrices allow us to build the coupling ma¬ 
trix V, the self-energies E, and the surface Green’s functions 

@s, o : 

V = K-! (T _1 - T) 

E = R\ T 

G*s,o = -(^o + S)- 1 (20) 


FIG. 4: (Color online) Infinite system whose unit cell is the EM re¬ 
gion. This is linked to neighboring EM cells by periodic boundary 
conditions. 
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FIG. 5: (Color online) Super cell containing the extended scattering 
region EM. The EM region is surrounded by buffer vacuum regions 
to its left and right. 


The procedure described above to compute the lead Green’s 
functions can fail, because of the singular behavior of the 
Hamiltonians matrices K\, which lead to numerical inaccu¬ 
racies in the solution of Eq. and is usually manifested 
in the program producing a number of positive and negative 
channels different from N. Notice that if the number of pos¬ 
itive and negative channels is different from N, then the dual 
vectors cannot be found by inverting Q and Q. There exist 
several schemes to regularize K\, based on decimating out 
the offending degrees of freedom. These procedures are ex¬ 
plained in detail in Refs. 1 91371 1. GOLLUM uses a suitable 
adaptation of these methods, which is described in the ap¬ 
pendix. 


3. Generating the Hamiltonian of the extended scattering region 

j^EM 

As noted above, JC EM can be provided as a model Hamil¬ 
tonian, or generated by a DFT or other material-specific pro¬ 
gram. One of the strengths of GOLLUM is an ability to treat 
interfaces with high accuracy. In a tight-binding description, 
tight-binding parameters of a particular material are often cho¬ 
sen by fitting to a band structure. However this does not solve 
the problem of choosing parameters to describe the interface 
between two materials. Often this problem is finessed by 
choosing interface parameters to be a combination of pure- 
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material parameters such as an arithmetic or geometric mean, 
but there is no fundamental justification for such approxima¬ 
tions. 

Therefore we describe here methods to generate JC EM us¬ 
ing a DFT program, where the inclusion of branches as part 
of the extended scatterer occurs naturally. Fig. 0 shows 
an example of a junction where the electrodes are identical. 
The system is composed of super cells formed from a central 
scatterer and PLs. There are periodic boundary conditions in 
the longitudinal direction, such that the TPL of one branch of 
a super-cell is linked smoothly to the TPL of a neighboring 
super-cell. Running a DFT program for such a super-cell then 
automatically generates JC EM . Provided the super cells con¬ 
tain sufficient PLs, the Hamiltonians K PPL and K PPL as¬ 
sociated with the TPLs will be almost identical to those gen¬ 
erated from a calculation involving an infinite periodic lead, 
ie. if the Hamiltonians Kq and K\ associated with the PLs 
are generated from a calculation involving an infinite peri¬ 
odic lead, then provided the super cells contain sufficient PLs, 
these will be almost identical to K PPL and K PPL respec¬ 
tively. In this case then there will be minimal scattering caused 
by the junction between the TPL and the lead. Clearly there 
is a trade-off between accuracy and CPU time, because insert¬ 
ing more PLs increases the size and cost of the calculation. In 
practice, the number of PLs retained in such a super-cell is in¬ 
creased in stages until the results do not change significantly 
as the number of PLs is increased further. 

There exist situations where the electrodes are dissimi¬ 
lar, either chemically, or because of their different crystalline 
structure, or because their magnetic moments are not aligned. 
In these cases, there cannot be a smooth matching between 
TPLs of neighboring super cells in Fig. 4. To address this 
situation, we use a setup similar to that displayed in Fig. 0, 
where additional PLs are appended to the branches in the EM 
region. These additional PLs are terminated by artificial sur¬ 
faces and surrounded by vacuum. The TPLs are then cho¬ 
sen to be one of the PLs near the middle of each branch and 
should be surrounded by enough PLs both towards its artifi¬ 
cial surface and towards the central scattering region. Then 
the PLs placed between the TPL and the artificial surface are 
discarded. These sacrificial PLs ensure that the chosen TPL is 
unaffected by the presence of the artificial vacuum boundary. 
Clearly, calculations of this sort are more expensive in numer¬ 
ical terms than those performed with super cells generated as 
in Fig. 0, because they contain many more atoms. 


4. Hamiltonian assembly 


In an ab initio calculation of the transport properties of a 
junction, the DFT program produces the Hamiltonians and 
Fermi energy of the EM region 'H EM , S EM and E EM , and 
of each lead Hq ±1 , Sq ±1 and E F in separate runs. Notice 
that the Hartree potential is defined up to a constant, which is 
usually different for the EM and for each lead. This usually 
means that the energy origin of the EM and of the correspond¬ 
ing lead PLs, as well as their Fermi energies do not agree with 


each other, so Eq. (2) must be rewritten as follows: 

Kl ±l = Hl ±1 - (£*, + E) S 0 ,±i (21) 

where we have referred the energy of each lead to its own 
Fermi energy. To fix the Hamiltonian mismatch we define a 
realignment variable for each lead as follows: 

A* = /j,) — fj) (22) 

where /i indicates a relevant orbital or group of orbitals. Then, 
the Hamiltonian of each lead is realigned with that of the EM 

K'o,± 1 = K ±1 - (E™ + E + A*) 5 0 ,±1 

= flS,±i - (E l F + E) S 0 ,±1 (23) 

It turns out that the renormalized E F and bare E F Fermi en¬ 
ergies of each lead do not match perfectly with each other if 
the number of PLs in the EM region is not sufficiently large. 
This is the case when for efficiency reasons, it is desirable to 
artificially minimize the size of the EM. Sometimes it is ad¬ 
visable to choose the Fermi energy of one of the leads E F as 
the reference energy. In this case, a second overall shift can 
be performed using either A 1 or the quantity S = E F — E em . 


5. Scattering matrix and transmission in multi-terminal devices 


We note that the most general scattering state in a given lead 
i at a given energy E can be written as a linear combination 
of open and closed channels as follows 


\¥(E)) 


E ° k ’ 

hi 


IT' 


J Qi 


\EM. 


+ IV) (24) 


where ki and qi denote here open positive and negative chan¬ 
nels and | ^ l (ki)) and 1are their normalized kets. Here, 
the contribution of all the closed channels in lead i is described 
by the ket \x l ) • Consequently, the number of electrons per unit 
time flowing between two adjacent PLs within the lead is 


r(E)= y> fe j 2 - E>i 2 (25) 

hi£i qi£i 


We pick in this section the convention that positive direction 
in the lead means flow towards the EM region and vice versa. 
So positive (negative) open channels are also called incom¬ 
ing (outgoing) channels of lead i. With this notation, the 
wave-function coefficients of the incoming open channels of 
a given lead are determined by the properties of the reservoir 
connected to the lead. 

The wave function coefficients of the open outgoing chan¬ 
nels of lead i are obtained from the amplitudes of all incoming 
channels by 


4 = E 44 4 (26) 

j,k 

where qi ( kj ) is an outgoing (incoming) dimensionless wave- 
vector of lead i (j). It is therefore convenient to assemble the 
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wave-functions of the M l outgoing and M l incoming open 
channels of a given lead i in the column vectors (5\ 0\ and 
all the scattering matrix elements connecting leads i and j into 
the matrix block S l E Notice that the dimensions of S l i are 
x . Then the above equation can be written more com¬ 
pactly as 


o' \ 


( S 11 S 12 ... S 1P \ 


( o 1 ^ 

o 2 

= 

g 21 g 22 s 2 P 


o 2 

U p J 


{ S P1 S P2 ... S pp ) 


\o p ) 


By normalizing the Bloch eigenvectors C(k), C(q ) and 
their duals to unit flux, 


where 


KU = - 


E l Jo 

1 Qcxikp. 


(33) 


and (E) is the number of open incoming channels of type 
a , energy E in lead i. Note that in the above summation, 
q ai runs over all outgoing wave-vectors of energy E and type 
di of lead i and kp j runs over all incoming wave-vectors of 
energy E and type fa in lead j. 

If i and j are different leads, then sq.^ is often called the 
transmission amplitude and denoted tq ^., while if they are the 
same lead, then Sq^. is called the reflection amplitude 
Similarly, for i ^ j, it is common to define the transmission 
coefficient ^ as 


C\ki) = C^fa)/^ V\ki) = D\ki) 

c'fe) = cr^/y/v^, z>‘(©) = Di tii) (28) 


rpij 
1 &i,dj 


x ; i4,*„i 2 

Qati 


(34) 


the matrix elements of the scattering matrix block connecting 
leads i and j can be written as. 

4k t = &(®) ( V 4 - I Sij ) C? fo) (29) 


Here G 3 g is the off-diagonal block of the surface Green’s func¬ 
tion defined in Eqs. (11) and (12), that connects leads i and j 
and V 1 is the matrix defined in Eq. (20). 

With the above notation, if the incoming channel fa of lead 
i is occupied with probability /&. (. E ) (ie if in Eq. (25), o/ Ci = 
1 with probability fa .(E)) then the number of electrons per 
unit time, entering the scattering region from reservoir i along 
channel ki with energy between E and E + dE is 


and for i = j, we define the reflection coefficient as 


R 


oti,di 


£ i 



(35) 


so that 

dE 

dlU E ) = -r- { £ [ ^* (£;)5a * ft ^‘ R “i.A ]/ A( £7) 

di 

- £ Efafim ( 36 ) 

j^hdj 

Note that unitarity of the scattering matrix requires 


dI%(E) = (dE/h)f ki (E) (30) 

and the number per unit time, per unit energy leaving the scat- 
terer and entering reservoir i along channel q h with energy 
between E and E + dE is 


. E !*£,*./= E i*SU/= 1 (37) 

iQ-cti Oidjkpj 

Hence the sum of the elements of each row and column of 
the matrix P is zero: 


dI^(E) = (dE/h) J 2 1 4 kffkA E ) (3D 

QiJ,kj 


E pij = pi,3 = n 

^<Xi,dj Z . U 

j,/3j i, a i 


(38) 


In many cases, the incoming and outgoing channels of each 
lead i can be grouped into channels possessing particular at¬ 
tributes (ie quantum numbers) labeled c^, fa .etc. This 

occurs when all incoming channels of a particular type oli in 
lead i possess the same occupation probability f&(E). For 
example, all quasi-particles of type in reservoir i may pos¬ 
sess a common chemical potential p (Xi and ffa (E) may take 
the form ffa(E) = f(E — /x a .), where f(E) is the Fermi 
function. In this case, if the incoming and outgoing channels 
of type oti belonging to lead i possess wave-vectors k a . 9 q a . 9 
then the number of quasi-particles per unit time of type oti 
leaving reservoir i with energy between E and E + dE is 

dP ai (E) = ( dE/h ) £ 4 (E) (32) 

Old j 


or equivalently, 


£^U + £ Efa < 

di j^hdj 


and 


E*& 


/ -J CX i 5 dj 




(39) 


(40) 


From Eqs. (36) and (39), if fp.(E) is independent of j 
and fa then dlfa(E) = 0 for all i and as expected. For 
this reason, in the above equations, ( E ) can be replaced by 
f J p. (E) = (E)—f(E) 9 where /(E) is an arbitrary function 

of energy, which in practice is usually chosen to be a Fermi 
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function, evaluated at a convenient reference temperature and 
chemical potential. 

When comparing theory with experiment, we are usually 
interested in computing the flux of some quantity Q from a 
particular reservoir. From Eq. ( [32] ), if the amount of Q carried 
by quasi-particles of type is Q ai (E), then the flux of Q 
from reservoir i is 

r Q = f(dE/h ) £ Qai pp, (E) (41) 

In the simplest case of a normal conductor, choosing Q a . = 
—e, independent of oti , the above equation yields the electrical 
current from lead i. Within GOLLUM a;* may represent spin 
and in the presence of superconductivity may represent hole 
(ai = h) or particle (c^ = p ) degrees of freedom. In the 
latter case, the charge Q p carried by particles is -e, whereas 
the charge Qh carried by holes is +e. 

B. Incorporation of non-equilibrium effects in the 
transmission coefficients 

GOLLUM starts from a mean-field Hamiltonian provided 
either by the user or by an outside material-specific DFT code. 


It then computes the scattering matrix and its related trans¬ 
port properties. When finite voltages are applied to the elec¬ 
trodes, they change the distribution of incoming and outgo¬ 
ing electrons and therefore the underlying Hamiltonian. For 
example, a finite voltage in a two-terminal device may intro¬ 
duce an electrostatic potential, which should be included in 
the Hamiltonian. A key feature of many NEGF codes in¬ 
cluding SMEAGOL is that such effects can be treated self- 
consistently, albeit at the cost of a greatly increased computing 
overhead. To avoid this overhead, GOLLUM assumes that the 
user is able to provide a modified Hamiltonian at finite volt¬ 
ages. 

Based on our experience on the development and usage of 
NEGF programs and as demonstrated in in section III.J be¬ 
low, we have found that in many cases, the following intuitive 
modification of the initial zero-voltage Hamiltonian yields 
reasonable-accurate voltage-dependent transmission coeffi¬ 
cients Tij(E , V) connecting leads i and j. The scheme en¬ 
ables the simulation of non-trivial I — V curves which com¬ 
pare favorably to those obtained using NEGF techniques and 
enables the modeling of generic non-equilibrium transport 
phenomena such as negative differential resistance (NDR)jmd 
current rectification in close agreement with NEGF codeP^. 

Consider the case where each lead has a different voltage 
V 1 . Then the finite-voltage Hamiltonian takes the form 


JC = 


1C 4 


K M4 


eU 4 

s 4 


0 

0 



0 


/C 4M 

— eV 4 

£4M 

0 


K? - 

- eV 3 S 3 

0 



0 


/c 3M 

-eV 3 

<s 3M 

0 



0 

K 2 - eU 2 

s 2 


0 


AC 2M 

— eV 2 

<S 2M 

0 



0 

0 


K 1 - 

-eV 1 

s 1 

K 1M 

— eV 1 

<S 1M 

eV 4 

£M4 

1C M 3- 

- eV 3 S M3 

JC M2 - eU 2 

s M2 

/C M1 - 

-eV 1 

s M1 


AC em 



\ 


/ 


(42) 


We find that JC EM needs only be computed at zero voltage 
in most cases; the effect of a finite bias can be accounted for 
by a suitable re-alignment of the energies of the orbitals in 
the EM region with the shifted energy levels of the electrodes. 
Mathematically, we apply a simple shift to the Hamiltonian 
matrix elements at each orbital n in the EM region 

/C EM —> /C EM (V) = JC EM - eV n <S EM (43) 

where these local shifts V n depend on the junction electro¬ 
statics, which in many cases are known. For example, in the 
case of a highly-transparent junction, the shifts can be mod¬ 
eled by a linear voltage ramp connecting the matrix elements 
of the orbitals at the TPLs of the EM region. In contrast, when 
the central scattering region Hamiltonian is connected to 
the PL Hamiltonians Kq in each branch of the EM region 
by weaker links K", the voltage drop and therefore the re¬ 
sistance is concentrated at these spots. In this case, we take 
V n = y i for all orbitals in branch i , starting at the TPL and 
up to the linker atoms, and Vi m 0 for all the orbitals inside 


the M region itself. This scheme performs specially well for 
systems where the states around the Fermi level (HOMO or 
LUMO) are localized at or close to the contact atoms. It en¬ 
ables us to mimic accurately junctions displaying non-trivial 
negative differential resistance, as well as rectification effects 
for asymmetric molecules 38 . 


C. Virtual leads versus physical leads. 

What is the difference between a lead and a channel? From 
a mathematical viewpoint, channels connect an extended scat¬ 
tering region to a reservoir and the role of lead i is simply to 
label those channels which connect to a particular reser¬ 
voir i. Conceptually, this means that from the point of view of 
solving a scattering problem at energy E , a single lead with 
N(E) incoming channels can be regarded as N(E) virtual 
leads, each with a single channel. GOLLUM takes advantage 
of this equivalence by regarding the above groups of channels 
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with wave-vectors k ai ,q ai as virtual leads and treating them 
on the same footing as physical leads. From this viewpoint, 
Eq. ( [32] ) and ( [36] ) yield the number of quasi-particles per unit 
time ’’from virtual lead oti ” entering the scattering region with 
energy between E and E + dE. 

This viewpoint is particularly useful when the Hamiltonians 
Hq, H\ describing the PLs of the physical lead i are block di¬ 
agonal with respect to the quantum numbers associated with 
k ai , q ai . For example, this occurs when the leads possess a 
uniform magnetization, in which case the lead Hamiltonian is 
block diagonal with respect to the local magnetization axis of 
the lead and a represents the spin degree of freedom a. This 
occurs also when the leads are normal metals, but the scat¬ 
tering region contains one or more superconductors, in which 
case the lead Hamiltonian is block diagonal with respect to 
particle and hole degrees of freedom and a represents either 
particles p or holes h. More generally, in the presence of both 
magnetism and superconductivity, or combinations of singlet 
and triplet superconductivity, a would represent combinations 
of spin and particles and holes degrees of freedom. 

In all of these cases, Hq, H[ are block diagonal and it is 
convenient to identify virtual leads with each block, be¬ 
cause GOLLUM will compute the channels k ai ,q ai belong¬ 
ing to each block in separate calculations and therefore guar¬ 
antees that all such channels can be separately identified. This 
is advantageous, because if all channels of Hq, H\ were cal¬ 
culated simultaneously, then in the case of degeneracies, arbi¬ 
trary superpositions of channels with different quantum num¬ 
bers could result and therefore it would be necessary to imple¬ 
ment a separate unitary transformation to sort channels into 
the chosen quantum numbers. By treating each block as a vir¬ 
tual lead, this problem is avoided. Examples of this approach 
are presented below, when describing the scattering properties 
of magnetic or normal-superconducting-normal systems. 


D. Charge, spin and and thermal currents 


In the presence of non-collinear magnetic moments, pro¬ 
vided the lead Hamiltonians are block diagonal in spin indices 
(in general relative to lead-dependent magnetization axes) 


choosing ai = and Q ai = — e in Eq. (41) yields for the 
total electrical current 


Il = ~e 


J (dE/h) J2 P °t 


Jim 


(44) 




Note that in general it is necessary to retain the subscripts i, j 
associated with cy or oy, because the leads may possess dif¬ 
ferent magnetic axes. 

Similarly the thermal energy from reservoir i per unit time 
is 

p q = j (ds/h) ]r (E-umjjim (45) 


there is no spin-flip scattering, P^ 3 a , = P^ 3 a d a y . In this case, 
the total Hamiltonian of the whole system is block diagonal in 
spin indices and the scattering matrix can be obtained from 
separate calculations for each spin. We assume that initially 
the junction is in thermodynamic equilibrium, so that all reser¬ 
voirs possess the same chemical potential po. Subsequently, 
we apply to each reservoir i a different voltage Vi, so that its 
chemical potential is m = po — e Vi. Then from equation 
( [32] ), the charge per unit time per spin entering the scatterer 
from each lead can be written as 

K = (—e) [ (dE/h) PjfjE) (46) 

<J,j 

and the thermal energy per spin per unit time is 

I\ = f (dE/h) J2(E - hi)PjJi(E) (47) 

< 7,3 


where e = |e| and p a (E) = f(E - p, t ) - f(E - n) is the 

deviation in Fermi distribution of lead i from the reference 
distribution f(E — fi). 

In the limit of small potential differences or small differ¬ 
ences in reservoir temperatures, the deviations in the distri¬ 
butions from the reference distribution f 3 a (E) can be approx¬ 
imated by differentials and therefore to evaluate currents, in 
the presence of collinear magnetism, GOLLUM provides the 
following spin-dependent integrals 


L ijA T ) = /” dE (E — n 0 ) n (E, V = 0) 

(48) 

In the presence of two leads labeled i = 1,2, the spin- 
dependent low-voltage electrical conductance G(T), the ther¬ 
mopower (Seebeck coefficient) S e (T ), the Peltier coefficient 
n(T) and the thermal conductance k(T) can be obtained as 


G*(T) 

Sl(T) 

n AT) 

k<j(T) 


(e 2 /h)L° 12 , c 

1 Lja 


TSl(T) 

1 (r2 

hT [ L ^ L y ) 


(49) 


so that the equivalent spin-summed magnitudes are 


G(T) = ^ G a (T) 

<J 

S e (T) = J2 S t(T) 

<J 

n(T) = y^rUT) 

<j 

k(T) = y> CT (T) (50) 


For the special case of a normal multi-terminal junction 
having collinear magnetic moments, a% = a for all i and since 


Note that the thermal conductance is guaranteed to be posi¬ 
tive, because the expectation value of the square of a variable 





12 


is greater than or equal to the square of the expectation value. 
For a two-terminal system, the above expressions allow us to 
obtain the electronic contribution to the thermoelectric figure 
of meri@ 


ZT = 


1 


M 2 M 2 

(M 2 ) 2 


- 1 


(51) 


E. Additional functionalities 


1. Spectral adjustment 


A phenomenological scheme that improves the agreement 
between theoretical simulations and experiments in, for exam¬ 
ple, single-molecule electronics consists of shifting the occu¬ 
pied and unoccupied levels of the M region downwards and 
upwards respectively to increase the energy gajPEl of the M 
region. The procedure is conveniently called spectral adjust¬ 
ment in nanoscale transport (SAINT). At the request of a user, 
GOLLUM modifies the Hamiltonian operator of the M region 
as follows: 

k u = k° M + a q £1* no X* no | T~ A u ^ ^ |^nu X* nu | 

no nu 

(52) 

where A Q5ll are energy shifts and (no, nu) denote the occupied 
and unoccupied states, respectively. By using the definition of 
the density matrix operator of the M region 


PM = £ |*no)(^no| 
no 

t = £l*»X*n| (53) 

n 

The above Hamiltonian can be rewritten as 

K u = Km + (Ao - A u ) p M + A u I (54) 
The equation can also be written in matrix form as 

K m = K^ + (A 0 - A u ) Sm Pm Sm + A u Sm (55) 


To find the density matrix, we first solve the generalized 
eigenvalue problem: 


r — 
ri M L n ~ 

P ~ 

R = 
rfK^R = 



V&n Sm c „ 


(c' 1 ,...,c' n ) 
£°m ~EIp 


(56) 

(57) 

(58) 

(59) 


where Ip is the P x P identity matrix, and we have arranged 
the eigen-energies e° into a diagonal matrix e ° M . Then 


(Pm)/x,/z' — C n,/i C n,/i' (60) 


In the simplest case, for a single-molecule junction, the shifts 
A o, u are chosen to align the highest occupied and lowest un¬ 
occupied molecular orbitals (ie the HOMO and LUMO) with 
(minus) the ionization potential (IP) and electron affinity (EA) 
of the isolated molecule 


A° = cromo + IP 

A° = — (clumo+^A) (61) 


However the Coulomb interactions in the isolated molecule 
are screened if the molecule is placed in close proximity to 
the metallic electrodes. Currently, GOLLUM takes this effect 
by using a simple image charge model 40 , where the molecule 
is replaced by a point charge located at the middle point of the 
molecule and where the image planes are placed 1 A above the 
electrodes’ surfaces. Then the shifts are corrected by screen¬ 
ing effects as follows: 


A 0 


A 


u 


A° + 


A°- 


e 2 In 2 
8 ttcq a 
e 2 In 2 
Sir eo a 


(62) 


where a is the distance between the image plane and the point 
image charge. 


2. Coulomb blockade and Kondo physics 

Many nanoscale-scale junctions are expected to show 
Coulomb blockade behavior, and in specific situations also 
Kondo features 45 . These features can be demonstrated by gat¬ 
ing the junction, and should appear as Coulomb and Kondo 
diamond lines in contour density plots of the low-voltage 
conductance as a function of bias and gate voltages. These 
strong correlation effects are completely missing in conven¬ 
tional DFT. Accurate parametrizations of the ground-state en¬ 
ergy density functional of the single-channel Anderson model 
exist that allow a correct description of those phenomena 4 ® 7 . 
However, most nanojunctions are better modeled in terms of 
a multi-channel Anderson model as we have chosen to do in 
GOLLUM. This model is described by the Hamiltonian 


^And — ^Leads P Km P ^Coupling (63) 

where the Hamiltonian at the central scattering region is given 
by 


ILm ^ ^ ^In^m T - U ^ ^ ft m fli (64) 

m m>l 

Here the m -sum runs over the M correlated degrees of free¬ 
dom and includes the spin index, e m denote the on-site en¬ 
ergies and U is the electronic Coulomb repulsion, that is as¬ 
sumed to be the same for all degrees of freedom. 

We map the central scattering region Hamiltonian K ^ in 
Eq. (4) into KmoI to extract the self-energies Sm of the cor¬ 
related degrees of freedom, in the spirit of Dynamical Mean 
Field Theory™. Following Ref. (09]), our correlated degrees 


no 
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of freedom are a subset of the eigenstates of K^, that we will 
call here molecular orbitals. This contrast with the approach 
followed in Ref. f48| where the correlated degrees of freedom 
where taken to be atomic orbitals of transition metal atoms. 
For a four-lead device, the details of our implementation are 
as follows. We take the EM Hamiltonian 


/ ZfO 

T^mi 

7Tm2 

7Tm3 

7Tm4 ^ 

TGm 

Tfbranchl 

0 

0 

0 

7^2M 

0 

Tfbranch2 

0 

0 

^3M 

0 

0 

Tfbranch3 

0 

\ 7T 4M 

0 

0 

0 

Tfbranch4 J 


(65) 

where K ^ has dimensions P x P and therefore describes P 
molecular orbitals. We first solve the generalized eigenvalue 
problem at the M region as in Eqs. (54-57) in the previous 
section to find the rotation matrix R that diagonalizes : 

R'K&R = £°m-EI p (66) 


where I p is the P x P identity matrix. We then perform the 
direct product U = R <g> I to enlarge the size of the matrix 
R to the dimensions of the EM Hamiltonian matrix. We then 
rotate the EM Hamiltonian and compute the Green’s function 
of the EM region 

JC EM ' = UlC EMfi U ] (67) 



We now choose which molecular orbitals i = 1, ...M with as¬ 
sociated on-site energies e 0 are the correlated degrees of free¬ 
dom, and use the projectors Pi to find their projected Green’s 
functions and occupancies 

g\ = PiG EM ' Pi (69) 

V DFT = ~l I d£Imag[«?']/(£ F ) (70) 


where E F is the EM Fermi level. We then shift the on-site 
energy of the correlated orbitals by the conventional double 
counting term 50 51 


j^EM ' 


e l = e °-U(N EET -l/2) 


(71) 

(72) 


As a consequence, we have to recompute again the Green’s 
functions 


gEM 



9i = p i £ EM Pi 


(K co up ) f G^o AT coup ) 


(73) 

(74) 


which allows us to extract the hybridization function A°. 

The initial ingredients in the solution of the multichannel 
Anderson model are the set (e$, A°, U). They allow us to ex¬ 
tract the self-energies D*(e*, A°, U) using a impurity solver. 
These are added again to the on-site energies 

e* e* + £< (76) 

leading to a new EM Hamiltonian JC EM and associated 
Green’s function Q EM . From here we compute a new hy¬ 
bridization function 

A i = E-Si-- (77) 

9i 

with which new self-energies £* are determined. The cycle is 
repeated until self-consistency in A^ and is achieved. The 
resulting JC EM is inserted back into Eq. (12) and the surface 
Green-function matrix Gs is computed to extract the transport 
properties of the correlated junction. 

We have decided to include in GOLLUM a finite- V impu¬ 
rity solver. This way, we can subtract the double-counting 
terms and place the molecular orbitals at their correct bare 
energy positions by using Eq. CD- There exist a variety of 
finite -U impurity solvers based on perturbation expansions on 
the Coulomb interaction U, on the hybridization function A°, 
on interpolative approaches, on Monte-Carlo algorithms (see 
Ref. 08} for a detailed account of some of these solvers), 
or on Numerical Renormalization Group techniques 52 (NRG). 
NRG techniques have superior accuracy, but they bring high 
computational demands. Slave-boson-based expansions on 
A° like the OCApEH are rather accurate and less expensive 
numerically. 

The impurity solver used in GOLLUM is based on the Inter¬ 
polative Perturbation Theory approach!®^, where the second- 
order in U expression for the electron self-energy is interpo¬ 
lated to match the atomic self-energy, and adjusted to satisfy 
consistency equations for the high-energy moments together 
with Luttinger’s theorem. This approach is computationally 
very simple, but has been proven to provide reasonable re¬ 
sults for the multi-channel finite -U Anderson model® 7 . Its 
main shortcoming is that it overestimates the Kondo Temper¬ 
ature, as we discuss in Section Specifically, the im¬ 

purity solver that we have implemented to handle the multi¬ 
channel Hamiltonian ( [64] ) is described in Ref. (f55k although 
we have corrected errors in some of the equations in that refer¬ 
ence. We note that this impurity solver handles M > 2 spin- 
degenerate correlated degrees of freedom, so that M must be 
an even number. In other words, these channels must come as 
Kramers pairs. We stress that other impurity solvers can be 
implemented straightforwardly, due to the modular nature of 
GOLLUM. 


3. Inclusion of a Gauge field 


This is cast in the form 
1 

9i 


E-u- A? 


A? = E 


1 

9i 


To compute transport properties in the presence of a mag¬ 
netic field GOLLUM allows the user to introduce a Peierls 
substitution by changing the phase factors of the coupling 
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where Ii e f t (Iright) is the current from the left (right) reser¬ 
voir, fii e ft — /i (fright ~ AO is the difference between the 
chemical potential of the left (right) reservoir and the chem¬ 
ical potential fi of the superconducting condensate and the 
voltage difference between the left and right reservoirs is 
{fiieft ~ /iright)/ e. This expression is the low-voltage limit of 
more general current-voltage relations discussed in ll59l60l . 
The generalization to multi-probe structures is described in 
Refs. 1611621 , to thermoelectric properties of superconduct¬ 
ing nanostructures in Refs. 1631641 and to ferromagnetic¬ 
superconducting structures in Refs. I65H671. In this equation, 


FIG. 6: Two-probe device consist of reservoirs a and [3 connected to 
a superconductor 


/ Mi e f t - R 0 + R a —T' 0 + T' a \ 
|| — T 0 + T a M ri gh t — R! 0 + R' a J 


elements 58 between atomic orbitals. For example in the case 
of a nearest-neighbor tight-binding Hamiltonian, the inter- site 
matrix element Hij between site i and site j is replaced with 
the modified element, 


K = H io* 


-Hfr r jMr)dr 


(78) 


where and r j are the positions of site i and j and A is the 
vector potential. The gauge is chosen such that the principal 
layers of the leads remain translational invariant after the sub¬ 
stitution. As an example, below we demonstrate how GOL- 
LUM describes the quantum Hall effect in a disordered square 
lattice, with a perpendicular uniform magnetic field. 


4. Superconducting systems 

Figure |6fa) shows a two-probe normal-superconductor- 
normal (N-S-N) device with left and right normal reservoirs 
connected to a scattering region containing one or more super¬ 
conductors. If the complete Hamiltonian describing a normal 
system of the type shown in Fig. 2 is Hjg, then in the presence 
of superconductivity within the extended scattering region, the 
new system is described by the Bogoliubov-de Gennes Hamil¬ 
tonian 



where the elements of the matrix A are non-zero only in the 
region occupied by a superconductor, as indicated in Figure 
[6jb). Physically, Hn describes particle degrees of freedom, 
describes hole degrees of freedom and A is the super¬ 
conducting order parameter. 

The multi-channel scattering theory for such a normal¬ 
superconducting-normal (N-S-N) structure was first derived 
by Lambert in Ref. (59|, where the following current-voltage 
relation was presented: 

| Iieft | 2 e 2 / {/iieft /i) /& 

y Iright J ^ \ {/bright a0/^ 



where Mi e f t ( M rig ht) is the number of open channels in the 
left (right) lead, R 0 , T 0 (. R a ,T a ) are normal (Andreev) reflec¬ 
tion and transmission coefficients for quasi-particles emitted 
from the right lead, R' 0 , T' 0 ( R ' a , T' a ) are normal (Andreev) re¬ 
flection and transmission coefficients from the left lead and 
all quantities are evaluated at the Fermi energy E = pt. As a 
consequence of unitarity of the scattering matrix, these satisfy 


R 0 +T 0 +R a +T a = Mieft and R' 0 +T/-\-R'a+T^ = M r i g ht • 

The current-voltage relation of Equ. (80) is fundamen¬ 
tally different from that encountered for normal systems, be¬ 
cause unitarity of the s-matrix does not imply that the sum of 
each row or column of the matrix a is zero. Consequently, 
the currents do not automatically depend solely of the ap¬ 
plied voltage difference {/iieft ~ /iright) /e (or more gener¬ 
ally on the differences between incoming quasi-article distri¬ 
butions). In practice such a dependence arises only after the 
chemical potential of the superconductor adjusts itself self- 
consistently to ensure that the current from the left reservoir 
is equal to the current entering the right reservoir. Insisting 
that I le ft = -Iright = L then yields 


2 e 2 / {/iieft fi) /& 

^ {/iright AO/^ 


= a 


-l 


I 

-/ 


(82) 


and therefore the two-probe conductance G = I/{{in e f t — 
fright) /e) takes the form of 


G = 


2 e 2 


^ 11^22 — ^ 12^21 


h an H - ^22 H - &12 + a 2i 


(83) 


The above equation demonstrates why a superconductor 
possesses zero resistivity, because if the superconductor is dis¬ 
ordered, then as the length L of the superconductor increases, 
all transmission coefficients will vanish. In this limit, the 
above equation reduces to {h/2e 2 )G = 2/R a + 2/R' a . In con¬ 
trast with a normal scatterer, this shows that in the presence of 
Andreev scattering, as L tends to infinity, the resistance ( = 
1/conductance) remains finite and therefore the resistivity (ie 
resistance per unit length) vanishes. 

In the notation of Eqs. ([32]) and ([36]), the above current- 


voltage relations and their finite-temperature, finite voltage 
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generalizations can be obtained from Eq. ( [41] ) by writing 
di = pi or hi to yield 


poo 

i l e = J ( dE/h ) E E E qmkUAw 

0 ai=pi,hi j=l,2 /3 j =p j ,h j 

(84) 

Since Q Pi = — e and Q\ li = +e, this becomes 




/•OO 

n = (-e/h) dE e E 

j «t, 2 (3j =Pj , hj 

(85) 

Since, in the low-bias limit, /^.(E) = —//.(E) = 
(—df(E)/dE)(/jbj — fi), where /(E) is the Fermi distribution 
with chemical potential /i, this simplifies to 


(e 2 /*) E “ M) 

j=l,2 


( 86 ) 


where 


— 


POO 

/ dE(-df(E)/dE) E KU ( 87 > 


ft3 ~Pj 5 h j 


The total current is obtained by multiplying Eq. ( [85] ) by a 
factor of 2 to account for spin. On the other hand, in the limit 
of zero temperature, dE(—df(E)/dE) = 1/2 Hence in 
this limit, the current-voltage relation ( [85] ) reduces to Eq. (80). 


F. Multiscale tools 

Simulation of the transport properties of a nanoscale-scale 
junction involves three distinct tasks. First, model geome¬ 
tries must be generated. Secondly, the Hamiltonian for each 
geometry must be constructed. Thirdly, the s-matrix can be 
calculated and transport properties of the junction calculated. 
GOLLUM separates these three tasks into three different pro¬ 
cesses. An overview of the work-flow of a generic GOLLUM 
calculation is shown in Figure [7] The three consecutive stages 
of the work process are denoted by the three dotted rectan¬ 
gular boxes. The initial step consists usually of modeling 
the atomistic arrangement of the junction. An initial struc¬ 
ture is usually guessed, followed by geometry optimization 
or molecular dynamics simulations to obtain a more realistic 
atomic arrangement. This task can be performed by either ab- 
initio or classical molecular-dynamics methods. For systems 
containing a few hundred atoms, a quantum-mechanical DFT- 
based simulation is usually the method of choice. However, 
experiments are often performed under ambient conditions or 
in a liquid environment. In these cases that the microscopic 
model should include the atomic structure of the environment, 
as we show below in section III.GOLLUM addresses this task 
by using classical molecular dynamics to model the environ¬ 
ment and in the spirit of the Born-Oppenheimer approxima¬ 
tion, feeding snapshots of the associated electrostatic field into 
the the DFT-based mean-field Hamiltonian. 

A similar approach is used to model the evolution of 
mechanically-controlled break junctions upon stretching, 


where the atomistic arrangement of the junction evolves 
slowly in time. In this case, if the same experiment is repeated 
a number of times, the junction geometry will be slightly dif¬ 
ferent each time. Therefore, a proper statistical analysis of 
the junction geometries is mandatory and calculations of the 
associated distribution of transmission coefficients is required. 
The task of generating junction geometries is also better suited 
for classical molecular dynamics situations. GOLLUM also 
facilitates the use of combined DFT and classical molecular 
dynamics approaches to gain accurate, yet quicker simulation 
result^. A non-comprehensive set of software tools is listed 
in Figure [7] Once the atomic arrangements are generated, 
these are fed into the second stage, where the Hamiltonian ma¬ 
trix is generated. This stage is in practice independent of the 
previous geometry construction and can be run separately, tak¬ 
ing only the output geometries of the first stage. The junction 
Hamiltonian can be generated using a variety of tools, some 
of which are listed in box II in Figure [7] A popular approach 
is the use of DFT codes that are able to write the Hamiltonian 
in a tight-binding language. In this way, model tight-binding 
Hamiltonians can also be easily generated. Other approaches 
involve the use of Slater-Booster or semi-empirical methods. 
In addition, GOLLUM has the ability to modify suitably these 
Hamiltonian matrices as discussed above. For example, the 
Hamiltonian matrix can be modified to include scissor cor¬ 
rections, Coulomb-blockade physics, a gate or bias voltage, a 
magnetic phase factor or a superconducting order parameter. 
Finally, stage III is the actual quantum transport calculation. 
This takes the Hamiltonian matrix as an input and calculates 
the s-matrix and associated physical quantities, such the elec¬ 
trical or spin current, the conductance, or the thermopower. 


III. DEMONSTRATOR CALCULATIONS 

In this section, we present a diversity of calculations, which 
demonstrate the broad capabilities of GOLLUM. For simplic¬ 
ity, we begin with a set of calculations on model Hamilto¬ 
nians, which demonstrate that GOLLUM can easily handle 
tight-binding models for a range of physical systems. We then 
move on to more material-specific calculations, in which the 
Hamiltonian is obtained from DFT. These include examples 
exhibiting Kondo physics, Coulomb blockade and non-linear, 
finite-voltage effects. Next we present more computationally 
challenging calculations involving van der Waals interactions, 
environmental effects and series of geometries associated with 
break-junction measurements. Finally an example of a quan¬ 
tum pump is presented, which requires access to the phase of 
scattering amplitudes. We define the conductance quantum 
Go = 2e 2 //i, that will be used frequently below. 


A. Simple one-dimensional tight-binding two and four 
terminal device 

As a first example, we consider a simple one-dimensional 
tight-binding chain containing a single orbital per PL and a 
single impurity orbital in the EM region, as shown in fig[8ja). 
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FIG. 7: Typical GOLLUM work-flow with various optional software tools 



FIG. 8 : Simple tight-binding one-dimensional (a) two probe (b) four 
probe systems containing a single orbital per PL and a single impu¬ 
rity orbital at the EM region. The parameters of the tight-binding 
model are so — 0 , £i — 1 , 7 = 1 and a — 1.5 (taken in arbitrary 
units). 


FIG. 9: Transmission and number of open channels for the simple 
tight-binding one-dimensional chain shown in Fig. [ 8 ja ) as a function 
of the energy. Energies are referred to the Fermi energy Ep and are 
given in units of 7 . 


all identical. 


We take the following parameters, that are given in arbitrary 
units. Within the leads, the site energies are £0 = 0, and 
the nearest neighbor couplings are — 7 . The impurity has a 
site energy e 1 = 1 and is coupled to the leads by a hopping 
element —a. Results are shown for 7 = 1 and a = 1.5. The 
transmission coefficient for this chain is shown in figure [9] 

As a second example, we consider the four-probe structure 
of Fig. [8jb), that shares the same set of parameters as the two- 
probe model above. The various transmission coefficients for 
this structure are shown in figure [10| By symmetry, these are 


B. The quantum Hall effect 

As an example of a quantum transport calculation with a 
magnetic field, we demonstrate the quantum Hall effect within 
the simple tight-binding square lattice shown in the inset of 
Fig. |TT| The lattice constant is set to a = 1 A. The onsite 
energies of the perfect lattice are e = 3.35 eV, the hopping 
integrals at zero magnetic field are 7 = 1 eV and the Fermi 
energy is set at zero. The red area in the figure denotes a 
disordered portion of the lattice. In this disordered area, the 
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FIG. 10: Transmission and number of open channels for the four- 
probe device shown in Fig. |8jb) that has four one-dimensional chain 
leads, as a function of the energy. Energies are referred to the Fermi 
energy Ef and are given in units of 7. 


onsite energies are randomly varied as e' = e + £, where £ is 
a random number distributed with uniform probability in the 
range (—0.2,0.2), (—0.4,0.4) and (—1,1) eV (red, green and 
blue dashed curves, respectively). 

The transport direction is chosen to be the y axis (e.g.: from 
bottom to top) while the x axis goes along the horizontal di¬ 
rection. To demonstrate the quantum Hall effect we introduce 
a homogeneous magnetic field perpendicular to the square lat¬ 
tice, pointing out of the paper, which is expressed in units of 
Bq = 6.58 x 10 4 Tesla. With this setup the vector poten¬ 
tial is chosen so that the lead remains translationally invariant 
along the y direction. This means that we implement a Peierls 
substitution of the form 


lij = l e 


j B (Vi-VjK* 


( 88 ) 


where X{ and y L are the coordinates of the site i. With this 
modified Hamiltonian the conductance calculated by GOL- 
LUM is shown in Fig. <mj- This clearly shows the presence 
of quantum Hall plateaus, which are resilient to the presence 
of disorder. 


C. Superconductivity 

As an example of scattering in the presence of superconduc¬ 
tivity, we now compute the electrical conductance of the N-S- 
N structure shown in Fig. © which contains two supercon¬ 
ducting regions with order parameters Ai and A 2 = Aie z6> . 
Such a structure is known as an Andreev interferometer and 
was first analyzed in Refs. 11691701 . where it was predicted that 
the electrical conductance is a periodic function of the order- 
parameter phase difference 6 , with period 2i r. At that time, 
this effect was completely missing from the more traditional 
quasi-classical description of superconductivity. When the 
missing terms were restored, good agreement between quasi- 
classical theory and scattering theory was obtained 71 . 



1/B (1/Bo) 


FIG. 11: Conductance G in units of Go as a function of inverse mag¬ 
netic field with various level of disorder. The magnetic field unit is 
set to Bo = 6.58 x 10 4 Tesla. The inset shows the square lattice used 
for the calculation. The black area denotes a perfect square lattice. 
The red area denotes a disordered portion of the lattice, where the 
inter-site distances are slight modified from 1 Ato perturb the phase 
contribution. The onsite energy for the regular lattice is e = 3.35 eV, 
the coupling with zero magnetic field is 7 = 1 eV and the Fermi en¬ 
ergy is chosen as zero. In the disordered area the onsite energy is ran¬ 
domly varied as e' = e + £, where £ is a random number distributed 
with uniform probability in the range (—0.2, 0.2), (—0.4, 0.4) and 
(—1,1) eV (red, green and blue dashed curves, respectively). 



FIG. 12: Two-terminal device consisting of two physical leads con¬ 
nected to a scattering region containing two superconductors with 
order parameters Ai and A 2 . The left (right) physical lead consists 
of two virtual leads pi and hi ( P 2 and / 12 ) carrying particle and hole 
channels respectively. 


In the following calculation, the Hamiltonian Hn of Eq. 
( [79] ) is simply a nearest neighbor tight-binding Hamiltonian 
on a square lattice, with diagonal elements So = 0 and 
nearest-neighbor couplings with 7 = 1 (in arbitrary units). 
Within the regions occupied by superconductor j , (where 
j = 1 or 2) the top (particle) sites are coupled to the bot¬ 
tom (hole) sites by A j, with | Aj| = 0.1 given in units of 7 . 
For 6 = 0, Figure ( f]~3j ) shows the energy dependence of the 
Andreev refection coefficient R a and the normal and Andreev 
transmission coefficients T 0 and T a respectively. The green 
line in Figure (\3\ represents the number of open channels in 
electron (hole) conducting leads. As expected, the Andreev 
reflection coefficient is large for small energies and decreases 
for energies above |Ai|. Substituting the values of these co¬ 
efficients at E = 0 into Eq. ( [83] ) and evaluating them for all 
6 yields the conductance versus 6 plot shown in Figure [14] 
As expected, the conductance is an oscillatory function of the 
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FIG. 13: Transmission coefficients ( R 0 (dot dashed line line), R a 
(solid blue line), T 0 (dashed line) and the number of open channels 
in the left lead for the device shown in Fig. |T2] ), as a function of 
the energy, as a function of the energy. Energies are referred to the 
Fermi energy Ef and are given in units of 7. 



FIG. 14: Two-probe conductance G in units of Go, for a N-S-N 
structure shown in Fig. as a function of the phase difference 
6 between the two order parameters. 

order-parameter phase difference 0 with period 27r. 

D. Non-collinear magnetism 

In this section, we compute the electrical conductance of 
the structure shown in Figure which we again describe 
using a simple tight-binding model of the form 

U f H]\f + M z M x iMy 
yM x +iMy h n -m z 

The Hamiltonian Hn is simply a nearest neighbor Hamilto¬ 
nian on a square lattice, with diagonal elements So = 0 and 
nearest-neighbor couplings with 7 = 1 (in arbitrary units) 
and for simplicity we choose M z = 0 everywhere. The sys¬ 
tems consists of two magnetic islands with magnetic moments 
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FIG. 15: Two terminal device consisting of two physical leads con¬ 
nected to a scattering region containing two ferromagnetic islands 
with magnetic moments (Mj,M^,0) and (M%,My,0) The left 
(right) physical lead consists of two virtual leads ti and ( ^2 and 
4 , 2 ) carrying up-spin and down-spin channels respectively. 



FIG. 16: Two-probe conductance G in units of Go, for the structure 
shown in shown in Fig. ( p~5| ) as a function of the angle 0 between the 
two magnetic moments. 


(Ml, My, 0) and (M^,My, 0), connected to non-magnetic 
leads. Choosing the Fermi energy to be E F = 0, and eval¬ 
uating Eq. ( [44] ) at zero temperature, Figure ( fib) ) shows the 
resulting electrical conductance as a function of the angle 0 
between the two magnetic moments. As expected, the con¬ 
ductance is an oscillatory function of the magnetic angle 9. 

Having discussed model systems described by simple tight- 
binding Hamiltonians, we now turn to more material-specific 
descriptions based on DFT. We will use the program SIESTA 
in most of the calculations below, and will provide many of 
the simulation parameters to help people to reproduce our cal¬ 
culations. 



FIG. 17: Sketch of the EM setups used in the calculation of spin- 
resolved transport through nickel electrodes that corresponds to the 
EM unit cell shown in Fig. (5). The scattering region is hown in 
light blue, the first PL is shown in greyish blue, and the second PL is 
shown in dark blue. The second PL is followed by vacuum. 
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FIG. 18: (Color online) Spin-resolved transmission coefficients as a 
function of energy for nickel-chain junctions. The different curves 
correspond to simulation with different levels of accuracy in the 
k p erp summations: (a), (b), (c) and (d) correspond to 1, 4, 16 and 
64 A;-points, respectively. (1) and (2) correspond to parallel and anti¬ 
parallel configurations, respectively. 


E. Spin polarized transport and magnetoresistance in nickel 
chains 

GOLLUM can describe voltage-dependent spin-polarized 
transport in spin-active junctions, made from a variety of met¬ 
als, including iron 2829 , platinum or palladium 72 . To demon¬ 
strate this, we describe here the voltage-dependent spin¬ 
filtering and magneto-resistive behavior of a two-terminal 
junction where (001) fee nickel electrodes with parallel (P) or 
anti-parallel (AP) spin orientations are connected by a nickel 
atomic chain 973 75 . Notice that because we may have elec¬ 
trodes with AP spin orientations, we are forced to use EM 
setups such as those shown in Fig. (5). We sketch in Fig. 03 
the Scattering Region used in the present calculation. It com¬ 
prises a 6-atom-long nickel chain, the electrodes surfaces and 
the two branches. The electrodes surfaces contain 2/3 atomic 
layers with 4 atoms each. The left and right branches contain 
two PFs that have 4 atoms each. The second PF is followed 
by vacuum. We have checked that the transport results in this 
example are reasonably converged if we choose PF1 as the 
TPF, which means that PF2 is sacrificial. 

To find the junction Hamiltonian, we use the program 
SIESTA. We use the Generalized Gradient approximation 
(GGA) functional 76 and take the theoretical GGA lattice con¬ 
stant of 3.45 A for the PFs as well as inter-atomic distances 
of 2.27 A along the chain. We have employed a single-^ (SZ) 
basis set to span the valence states and a mesh cutoff of 400 
Ry to define the real-space grid where the density, potential 
and matrix elements are calculated. 

To understand the spin-polarized transport properties of 
the junction, we will analise below the spin-dependent trans¬ 
mission coefficients T a (E ), together with the spin-dependent 
charge currents. These are computed using the approximate 


FIG. 19: (Color online) Charge current of the junction shown in Fig. 
( fT7| ) in the (a) P and (b) AP spin configurations as a function of the 
bias voltage applied to the junction. The different curves correspond 
to different numbers of k± -points. 

expression 

e r eV/2 

dE T a (E,V m 0) (90) 

™ J-eV/2 

The above approximation is quantitatively accurate for small 
enough bias voltages (< 0.5 V) and also shows the expected 
qualitative behavior at larger voltages. We have found that the 
transmission coefficients and currents depend sensitively on 
the number of transverse k± -points taken along the plane per¬ 
pendicular to the transport direction. As we will show below, 
we need to use at least 16 k± -points to achieve convergence. 
In other words, a r-point calculation provides a poor estimate 
of the transport properties of these junctions. 

We plot T a (E) as a function of the energy referred to the 
Fermi energy of the Scattering Region for P and AP spin 
orientations in Fig. ( p~8] ). The upper panel of the figure 
shows that the transmission coefficients for the P configura¬ 
tion are strongly spin-polarized. This polarization remains at 
the Fermi level, which suggests that these junctions could act 
as spin filters. The bottom panel of the figure shows the trans¬ 
mission coefficients for the AP configuration. The fact that 
these are different from those of the P spin orientation hints 
that these junctions could show significant GMR ratios. To 
quantify these statements, we compute the charge current of 
the junction in the P and AP configurations Ip and Iap and 
plot them in Fig. The figure shows that indeed these 

junctions show magnetorresistive properties. The figure also 
demonstrates that the currents depend on the number of Ap¬ 
points. To further give quantitative estimates of the spin ac¬ 
tivity of the junctions, we define the spin polarization in the P 
arrangement and the GMR ratio as 

Pp = 7p,t - /p,4 (91) 

GMR(%) = — ~ — x 100. (92) 

Iap 
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FIG. 20: (Color online) (a) Spin polarization of the current in the 
P configuration, measured in fiA and (b) GMR ratio plotted as a 
function of the bias voltage applied to the junction shown in Fig. 
0 - The different curves correspond to different numbers of Ap¬ 
points. 



FIG. 21: (Color online) A junction showing a molecule made of three 
phthalocyanine units connected via butadiyne linkers placed on top 
of two graphene electrodes separated by a physical gap of length 
17.265 A. The molecule is placed 3.4 Aabove the graphene sheets. 


Where /p ?cr are the spin-dependent currents for the P orienta¬ 
tion. We show in Fig. ( [20] ) these two magnitudes as a func¬ 
tions of the bias voltage applied to the junction. We indeed 
find large spin signals for these devices. Furthermore, the 
figure demonstrates that at least 16 k± -points are needed to 
achieve converged results. 


F. Simulation of a graphene-based junction using a van der 
Waals Density Functional 

GOLLUM can profit from the improved chemical accu¬ 
racy delivered by the most advanced density functionals. As 
an example, we discuss here the transport properties of the 
junction shown in Fig. where a single phthalocyanine 

trimer molecule bridges two graphene electrodes separated by 
a physical gap of length 17.265 A. These graphene sheets 
are armchair-terminated and passivated by hydrogen atoms. 
Periodic boundary conditions are applied in the two direc¬ 
tions across the graphene plane. The phthalocyanine units 


FIG. 22: (Color online) Transmission curves as a function of energy, 
referred to the Fermi energy of the Scattering Region. The different 
curves correspond to different electrode separations. The gap length 
is changed by removing or adding carbon layers. 


are linked by butadiyne chains. The planar anchors couple 
to the graphene via interaction of the 7r-clouds and therefore 
an accurate description of the chemical bonding and transport 
properties can only be achieved by the use of a van der Waals 
density functional. We use here the implementation of Dion 
et al. in the SIESTA program 77 78 . We have computed the 
Hamiltonian and overlap matrix elements using a double-zeta 
basis set for all the elements in the simulation, together with 
a grid fineness of 200 Rydberg. By minimizing the energy, 
we find that the molecule sits at a height of 3.4 A above the 
sheets. 


We have studied the impact of the length of the electrode 
gap on the transport properties of the junction by attaching ad¬ 
ditional armchair layers to the edges of both sheets; these lay¬ 
ers are made of two carbon rows, and have a width of 2.502 A. 
The transmission coefficients for several gap widths are shown 
in Fig. (22). We find several Breit-Wigner resonances associ¬ 
ated with molecular levels of the trimer. Remarkably, these do 
not shift in energy as the gap width varies 79 . However, deep 
dips appear for several gap widths. These are associated with 
interference among the different paths whereby electrons can 
propagate between the molecule and the electrodes. 


We have also studied the change in the transmission curves 
as the molecule is displaced laterally and longitudinally across 
the physical gap. We show representative examples of the 
transmission curves for longitudinal displacements in Fig. 
( [23]) . The figures show that the energy positions of the molec¬ 
ular Breit-Wigner resonances remain almost constant. We 
have found the same behavior for other graphene-based junc¬ 
tions: the energy position of the Breit-Wigner resonances 
for a given graphene-based junction does not depend on the 
molecule position relative to the physical gap, provided that 
the bonding mechanism is by physisorption. This universal¬ 
ity arises because physisorption carries no charge transfer be- 
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FIG. 23: (Color online) Transmission coefficient as a function of 
energy referred to the Fermi energy of the Scattering region. The 
different curves correspond to different longitudinal displacements 
of the molecule referred to the position shown in Fig. ([22]). The 
physical gap width is 14.763 A. 


tween the molecule and the sheets. Furthermore the elec¬ 
trodes are made from the same material and therefore there 
is no dipole moment associated with the contacts. Finally, 
the 7r — 7r hybridization between molecular orbitals and the 
electrode states is weaker than for the bonds present in most 
noble-metal/single-atom contacts, and does not have a large 
impact on the nature of the molecular orbitals. 


G. LDA+U description of gold porphyrin junctions 


As an example of a GQLLUM calculation using a LDA+U 
density functional® 81 } we describe in this section a case 
where strong electronic correlations may affect the trans¬ 
port properties of a nanoscale-scale junction. We discuss 


the junction shown in Fig. (24). Here gold (001) elec¬ 
trodes bridge either a porphyrin (P) or a metallo-porphyrin 
molecule (CuP or CoP). Electron flow through any of these 
three porphyrin molecules is carried by molecular orbitals that 
hybridize strongly with the gold s-orbitals. This gives rise 
to broad Breit-Wigner resonances in the transmission coef¬ 
ficients that are identical for the three molecules. However, 
for the CuP and CoP junctions, additional electron paths are 
created whereby electrons hop into and off the localized d- 
orbitals of the transition metal atom. The interference between 
direct and d-orbital-mediated paths creates sharp Fano reso¬ 
nances that can however be masked by the much wider Breit- 
Wigner resonances 82 . We see below how including strong cor¬ 
relations in the d-orbitals of the Co and Cu atoms in terms of 
a LDA+U approach produces strong shifts in the energy de¬ 
pendence of those resonances. 

We have computed the Hamiltonian using the SIESTA 
code, where we have picked a single zeta basis for the gold 
atoms at the electrodes, a double-zeta-polarized basis set for 
all the atoms in the molecule and a GGA functional. We have 
included a U correction term for the d-orbitals of the Cu and 
Co atoms in a mean-field fashion, in the spirit of the LDA+U 
approach 8081 . We present here our results for values of U 



FIG. 24: (Color online) Schematic view of a metallo-porphyrin 
molecule sandwiched by gold leads. Yellow, red, cyan, green, blue 
and orange represent gold, sulfur, hydrogen, carbon, nitrogen and Co 
or Cu atoms, respectively. 
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FIG. 25: (Color online) Transmission coefficient T(E) of the junc¬ 
tion shown in Fig. ( [24] ), where the central molecule is CuP. T is 
plotted as a function of the energy E referred to the Fermi energy 
Ef of the Scattering Region. The different panels correspond to the 
different U corrections added to the DFT Hamiltonian (see text). The 
green and gold ellipses encircle masked Fano resonances. They are 
originated by paths hopping onto the Cu 3d xz or d yz orbitals. The 
red ellipse circles a sharp Breit-Wigner resonance coming from C 
and N atoms. 


equal to 0, 2.5 and 4.5 eV. 

We find that the transmission coefficients of the three 
molecules display the same wide Breit-Wigner resonances, 
that correspond to molecular orbitals hybridizing strongly 
with the electrodes. These are shown in Figs. ( [25) and ( |26) 
for CuP and CoP respectively. In addition, the three molecules 
show a sharp Breit-Wigner resonance that is marked by a red 
ellipse in the figures. This resonance corresponds to a molecu¬ 
lar orbital encompassing C and N atoms that is weakly bonded 
to the electrodes. Interestingly, this sharp Breit-Wigner reso¬ 
nance shifts in energy if we change the value of the Coulomb 
interaction U for the CuP and CoP junctions. To understand 
this phenomenon, we have looked at the density of states of 
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FIG. 26: (Color online) Transmission coefficient T(E) of the junc¬ 
tion shown in Fig. \2A\ , where the central molecule is CoP. T is 
plotted as a function of the energy E referred to the Fermi energy 
Ff of the Scattering Region. The different panels correspond to 
the different U corrections added to the DFT Hamiltonian (see text). 
The green, gold and blue ellipses encircle Fano resonances. They 
are originated by paths whereby electrons hop onto the Co 3d x z or¬ 
bital (green ellipses) and from a Molecular Orbital composed by the 
Co 3d x y and 3d y z orbitals (blue ellipse). The red ellipse encircles a 
sharp Breit-Wigner resonance coming from C and N atoms. 


the junction projected onto each atomic orbital, and the local 
density of states integrated in a narrow energy window around 
the red resonance. We have found that the N- and C-based 
molecular orbital hybridizes with the d xy orbital of the cop¬ 
per or cobalt atoms and is therefore affected by the U -term. 
We have found additional sharp peaks appearing in T(E ) for 
the CuP and CoP junction that do not show up for the simple 
porphyrin junction. These are marked by green, blue and gold 
circles in Figs. 


and (26). These seem to be sharp Breit- 
Wigner resonances also. However, we have demonstrated 82 
that they are actually Fano resonances, where the Fano dip 
is masked by the transmission of the neighboring wide Breit- 
Wigner resonances. By plotting the density of states of the 
junction projected in each orbital, we indeed find that they cor¬ 
respond to the copper d xz or d yz orbitals. Because these Fano 
resonances are associated with atomic d-orbitals strongly lo¬ 
calized in the transition metal atom, we expect that adding 
a U -term will have a strong impact on their energy position. 
Fig. ( [25] ) shows how these resonances indeed shift in energy 
as U is increased. Note that one of the Fano resonances com¬ 
ing from the d xz copper orbital is strongly pinned to the Fermi 
energy, while other resonances rapidly move down in energy. 


H. Fixing the transport properties of OPE molecular junctions 
via the SAINT method 

We analyze in this section the transport properties of a se¬ 
ries (111) gold junctions that are bridged by OPE derivatives. 
The backbone of these molecules has a varying number of 
rings ranging from one to three. The molecules may be ori¬ 



FIG. 27: (Color online) A (111) gold junction sandwiching a tricene- 
dithiol molecule. The molecule is coupled to one ad-atom on one side 
and to a hollow gold site at the other end, and its orientation is titled 
with respect to the electrodes’ normal line. 


ented fully perpendicular to the electrodes surfaces, or making 
a tilting angle, as we show in Fig. ( [27] ) It is well established 
by n ov J — ^ that gold junctions that contain conjugated thiol- 
terminated molecules like OPEs have a larger conductance 
when the molecule is tilted. This is due to the increased over¬ 
lap of the p z states of the sulfur atoms when the angle between 
the molecule and the normal to the surface increases. We show 
in this section that a plain DFT-based calculation predicts that 
the largest conductance occurs when the molecule is oriented 
perpendicular to the electrodes. This deficiency is remedied 
by the use of the SAINT method. This method is an efficient 
semi-empirical correction that allows us to obtain quantitative 
agreement between DFT calculations and experiment^®®, 
as we have already stressed in Section II. 

The Hamiltonian of the junction has been obtained with 
the code SIESTA and a Focal Density approximation (FDA) 
functional. We have picked a single-zeta basis for the gold 
atoms of the electrodes, and a double-zeta-polarized basis for 
the atoms in the molecule. The PFs of the electrodes contain 
three atomic layers, each having 6x3 atoms. We have cho¬ 
sen junction geometries where the molecular derivatives are 
oriented either perpendicularly to the gold surfaces or making 
a 45 degrees angle, as shown in Fig. ( [27] ). Due to this tilt¬ 
ing angle, we had to use a non-periodic Scattering Region, as 
in Fig. (5). The scattering regions consisted therefore of the 
molecule, the two surfaces containing two atomic layers each 
and 3 PFs on each branch followed by vacuum. We chose PF2 
as the TPF, so all Hamiltonian matrix elements of PF3 were 
chopped off. 

We have computed the transmission curves T(E) for the 
referred OPE derivatives using conventional DFT, and have 
found for all of them that the conductance (computed from 
T(Ep)) is larger if the molecule is oriented perpendicular to 
the electrodes. The upper panel in Fig. ( [28] ) demonstrates 
this behavior for an OPE containing three rings. The figure 
shows that the higher conductance is originated by the posi¬ 
tion of the HOMO level of the molecule, that is placed only 
slightly below the Fermi energy of the Scattering Region. In 
contrast, the HOMO level of the tilted molecule is shifted far- 
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FIG. 28: (Color online) Transmission of the junction shown in Fig. 
{27}, containing an OPE molecule with 3 rings. T(E) has been cal¬ 
culated (a) without and (b) with SAINT corrections. Continuous 
and dashed lines correspond to perpendicular and tilted (45 degrees) 
molecules, respectively. 


FIG. 29: (Color online) Transmission of the junction shown in Fig. 

containing an OPE molecule with (a) 1, (b) 2 and (c) 3 rings. 
T(E) has been calculated with SAINT corrections. Continuous 
and dashed lines correspond to perpendicular and tilted (45 degrees) 
molecules, respectively. 


TABLE I: Corrections entering the SAINT scheme for occupied and 
unoccupied levels given in eV. 


# of Rings-tilt angle 

A 0 

A u 

1 -0° 

-1.7 

1.8 

1-45° 

-1.5 

1.6 

2 -0° 

-1.3 

1.3 

2-45° 

-1.2 

1.2 

U> 

O 

0 

-1.2 

1.3 

3-45° 

-1.1 

1.2 


ther away from Ep. This situation demands for the use of the 
SAINT correction scheme, that will reposition the molecular 
orbital levels at their correct energies. We show in the bot¬ 
tom panel of the figure that this is indeed the case, and that 
by the use of the SAINT scheme, the correct experimental 
trend is recovered, where tilted molecules show larger con¬ 
ductances. We have verified that the same change happens for 
OPE molecules containing one and two rings. Finally, we plot 
in Fig. (29) T(E) for the three molecules (containing one, 
two and three rings) for perpendicular and tilted orientations. 
( [30} we show the transmission of the OPE derivatives with a 
number of rings between 1 and 3 and two tilting angles, 0 and 
45 degrees. The figures show that all those junctions where 
the derivative is oriented perpendicularly to the gold surface 
(defined here to be 0 degrees) show a larger transmission at 
the Fermi level than the tilted cases in contrast with our ex¬ 
pectations discussed above 86 . 


The physical mechanism whereby the conductance of the 
tilted configuration is higher, is due to the higher hybridiza- 




FIG. 30: (Color online) Conductance of OPE molecules with a num¬ 
ber of rings between 1 and 3, calculated without (a) and with (b) 
SAINT corrections. Circles and squares correspond to perpendicular 
and tilted (45 degrees) molecules, respectively. 


tion between the molecular orbitals and the electrodes in the 
tilted configuration. This increases the width of the transmis¬ 
sion resonances and therefore decreases the effect of opening 
the gap with the SAINT correction. On the other hand, the im¬ 
age charge correction is also larger in the tilted configuration, 
since the molecule is closer to the surfaces, and therefore the 
reduction in the opening of the gap is also larger, which means 
the final gap ends up smaller in the tilted configuration. The 
conductance of each case is summarized in Fig. Notice 
that the SAINT correction scheme changes qualitatively the 
physical picture in this junction. 















































24 




FIG. 31: (Color online) Zero-voltage transmission coefficients T(E) 
of (a) a single-level symmetric Anderson model with input parame¬ 
ters defined in the main text; (b) a hydrogen atom sandwiched by 
(001) gold electrodes. The panels show curves computed at sev¬ 
eral temperatures to reflect the Kondo and the Coulomb blockade 
regimes. 


Our procedure for the SAINT correction scheme is as fol¬ 
lows. We first calculate the ionization potential (IP) and elec¬ 
tron affinity (EA) of the molecules in the gas phase. These gas 
phase corrections open the DFT HOMO-LUMO gap. How¬ 
ever, these bare shifts need to be corrected because of image 
charge effects. The final values for the correction shifts are 
summarized in table ^. Notice that the corrections are very 
similar in magnitude and have opposite signPl 


I. Kondo and Coulomb blockade effects 

As noted above, GOLLUM has a simple and flexible input 
data structure so that model Hamiltonians can be utilized eas¬ 
ily. As a simple example of a simulation exhibiting Kondo 
and Coulomb blockade behavior, we show here results ob¬ 
tained from GOLLUM for a tight-binding single-level An¬ 
derson model coupled to two semi-infinite chains, that cor¬ 
responds to taking M = 2 in Eq. ( [64] ). Due to the coupling to 
the two leads, the correlated level acquires a finite bandwidth 

V 2 

r = 27r V 2 p L = — (93) 

where pL is the density of states in the leads. For vanishing 
T the model is in the so-called atomic limit which is charac¬ 
terized by sharp peaks in the <i-level density of states pd at 
and 6 d + U. This limit corresponds to the Coulomb blockade 
regime in an actual junction where the conductance is strongly 
suppressed except at the charge degeneracy points. However, 
when the coupling to the leads increases (r becomes larger 
than the temperature T, but is still smaller than U), virtual 
processes allow the charge and spin in the molecule to fluctu¬ 
ate and a resonance close to the Fermi energy appears due to 
the Kondo effect. This simple model therefore captures rel¬ 
evant physics of molecular junctions such as the appearance 



FIG. 32: A hydrogen atom bridging (001) gold leads. The separation 
between H and each gold lead is 3.8 A. 


of the Coulomb blockade effect, and the crossover from the 
Coulomb blockade to the Kondo regime as the temperature is 
lowered below the Kondo temperature 




rF d (e d + U) I 
2 u r 


UT _ K u 

-e 8r 

2 


(94) 


where the last expression holds in the so-called symmetric 
limit €d + U/2 ~ 0. 

The interpolative impurity solver implemented in GOL¬ 
LUM provides a good quantitative description of the above 
phenomena in the weak coupling regime and is also qualita¬ 
tively correct in the intermediate and strong coupling regimes. 
It however overestimates the width of the Kondo resonance. 
We show in Fig. [31] (a) the zero-voltage transmission curve 
T(E ) computed with GOLLUM, and using the single-level 
Anderson Hamiltonian ( [64] ) in the symmetric limit. We take 
the following parameters: t *» 1 eV, U = 0.21 and V = 0.11 
so that T = 0.011 and 7rU/8r ^ 8, placing the junction in the 
strong correlation regime. These parameters yield a Kondo 
temperature Tk ~ 1.2 x 10 _5 t ~ 0.15 iT. The figure shows 
that the interpolative solution provides a transmission curve 
featuring the lower and upper Hubbard bands placed at their 
correct position and having the right width, together with a 
sharp Kondo resonance at low temperatures which progres¬ 
sively smoothens and eventually disappears as the tempera¬ 
ture is raised. However, the interpolative solution provides a 
Kondo temperature ~ 10 K, e.g.: two orders of magni¬ 
tude larger than the exact one. 

We now show the results obtained from GOLLUM for a 
similar junction, shown in Fig. ( [32] ), where a hydrogen atom 
bridges two gold (001) electrodes. In this case, the input 
Hamiltonian is generated by the DFT code SIESTA and the 
leads are repeated periodically in the plane perpendicular to 
the transport direction, using PL unit cells in AB stacking and 
3x3 atoms in each atomic layer. We have used a single-zeta 
basis set and a generalized gradient approximation functional. 
We have adjusted the distance d between the hydrogen atom 
and the leads to reproduce a coupling similar to that set for 
the above Anderson model which is achieved with d = 2.8 


A. The generated transmission curve is shown in Fig. 31 (b) 
for the same three temperatures used for the Anderson Hamil- 
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FIG. 33: Conductance as a function of gate voltage and temperature 
for a (001) gold junction bridged by a single hydrogen atom. 


tonian. We find that the shape of the transmission curves re¬ 
mains qualitatively the same. However, both the lower and up¬ 
per Hubbard bands and the Kondo resonance are now sharper 
at their tips. 

Using GOLLUM, we subsequently apply a gate voltage 
V g to the gold-hydrogen-gold junction and compute the low- 
voltage conductance G as a function of V g to compare the 
results using plain DFT versus DFT in combination with the 
interpolative method. The results are shown in Fig. (33). We 
have not included the double-counting term to make the com¬ 
parison between both approaches more explicit. The figure 
nicely shows how the interpolative method splits the single 
DFT peak into two Coulomb blockade peaks and also a Kondo 
peak. The figure also shows how the Kondo peak disappears at 
temperatures above leaving only the Coulomb blockage 
features. 

Finally, we subject the gold-hydrogen-gold junction to the 
combined effect of finite bias V and gate V g voltages. Fig. 
(34) shows density-contour plots of the low-voltage conduc¬ 
tance as a function of V and V g . This figure demonstrates 
that GOLLUM can nicely reproduce Coulomb blockade dia¬ 
monds, as well as the Kondo line, that disappears as the tem¬ 
perature is raised above . 


J. A junction displaying NDR behavior: a comparison 
between GOLLUM and SMEAGOL. 

We demonstrate with two examples how GOLLUM incor¬ 
porates finite-voltage effects. In the first, we have computed 
the current-voltage characteristics of a gold (001) junction 
sandwiching an alkane molecule, that we show in Fig. ( [35] ). 
We have computed the zero-bias Hamiltonian of the junction 
using the SIESTA code, with a double-zeta-polarized basis 
for all the atoms, and a GGA functional. The PLs contain 
two atomic layers, with 3x3 atoms each. We have applied 
periodic boundary conditions across the plane perpendicular 


to the transport direction and have computed H(k±) at the T 
point. 


For every given voltage V, we modify the EM Hamiltonian 
as described in Eq. All the orbitals n at the left branch 
in the EM region in Fig. (35) are shifted by V n = +U/2, 
starting at the TPL and stopping at the linking sulfur atom. 
Similarly, the orbitals at the right branch in the EM region are 
shifted by V n = — V/2 all, starting at the linking sulfur atom 
and including those at the TPL. The current I(V) is then com¬ 
puted from the modified Hamiltonian JC EM (V). The result¬ 
ing I — V curve is shown as a red dashed line in Fig. ( [36] ). 
For comparison, the black line in the same figure shows the 
I — V curves obtained with a full NEGF simulation using the 
code SMEAGOL. Finally, the dot-dashed green line shows the 
current-voltage curve obtained from GOLLUM by integrat¬ 
ing the zero-voltage transmissions obtained from JC EM (0). 
The figure demonstrates that our proposed method reproduces 
rather accurately the features found in the full NEGF calcu¬ 
lation, including the NDR feature at V « 2 volt, while the 
plain equilibrium calculation fails to reproduce the gross fea¬ 
tures of the current-voltage characteristics. It underestimates 
the low-voltage conductance by a factor of two. 


As a second example, we have calculated the current- 
voltage characteristics of a (111) gold junction sandwiching 
a porphyrin molecule. The junction geometry is similar to 
that shown in Fig. The sulfur atoms attach to the elec¬ 
trodes at a hollow site. The porphyrin molecule does not have 
in the present case a metallic atom at the center, but has two 
saturating hydrogen atoms instead. We have performed two 
series of calculations. In the first, the geometry and physi¬ 
cal gap distance has been relaxed and the calculations have 
been done at the most stable configuration. In the second, we 
have pulled the electrodes and sulfur atoms away. To do so, 
we have increased the sulfur-molecule distance by 0.3 A. The 
current-voltage curves are shown in Fig. These charac¬ 

teristics do not show non-trivial features. We note again that 
the plain calculation using zero-voltage transmissions fails to 
reproduce the NEGF curve, underestimating the conductance 
by a factor close to 2. In contrast, our prescription provides 
characteristics that reproduce accurately the NEGF results. 


To understand the difference between the dot-dashed green 
lines and the finite-voltage results of Figs. ([36]) and ([37]), we 


have analyzed the evolution of the transmission coefficients 
T(E,V) as the voltage bias V is ramped. We have found 
running NEGF simulations that the main non-equilibrium ef¬ 
fect that affects the current for the two junctions above is an 
energy shift of the molecular HOMO resonance, that moves 
up as the voltage is increased. Our prescription not only cap¬ 
tures the effect, but also follows accurately the evolution of 
the resonance shifts dictated by the NEGF calculation, at least 
at low voltages. By shifting the HOMO resonance upwards 
in energy, a larger weight of the resonance enters into the en¬ 
ergy integration window used to compute the current integral, 
hence increasing the current. In contrast, this effect can not be 
captured at all if one uses the plain T(E,V = 0) transmission 
coefficients. 
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FIG. 34: Density-contour plots of the low-voltage conductance G in units of Go of the gold-hydrogen-gold junction shown in Fig. ( [32] ). G is 
plotted as a function of the bias V in the vertical axis and the gate voltage V g in the horizontal axis. The conductance is plotted at a temperature 
(left panel) T — 0.11 K; (middle panel) T — 11 K; and (right panel) T — 1160 K. 



FIG. 35: (Color online) A junction where a Butane-dithiol mole 
is sandwiched by (001) gold electrodes and subjected to a finite 
potential. 



FIG. 36: Current-voltage curves of the junction shown in Fig. ( [35] ). 
The solid black line represents the result obtained from a full NEGF 
calculation using the code SMEAGOL. The red-dashed line shows 
the the curve obtained from GOLLUM using the method discussed 
in section II.B. The green dot-dashed line corresponds to integrating 
the zero-voltage transmission coefficient. 



y(V) 

FIG. 37: Current-voltage curves of a junction similar to that shown 
in Fig. ( [24] ), where gold (111) electrodes sandwich a porphyrin 
molecule, whose sulfur end-atom attaches to the electrodes at a hol¬ 
low position, (a) corresponds to the equilibrium distance between 
the electrodes and the molecule; (b) corresponds to a junction where 
each electrode and sulfur atom is pulled 0.3 angstrom away from the 
molecule backbone. The solid black line represents the result ob¬ 
tained from a full NEGF calculation using the code SMEAGOL. The 
red-dashed line shows the the curve obtained from GOLLUM using 
the method discussed in section II.B. The green dot-dashed line cor¬ 
responds to integrating the zero-voltage transmission coefficient. 


K. Temperature dependence of the thermoelectric properties 
of a C6o molecular junction. 


In this section, we show how GOLLUM can compute the 
thermoelectric properties of complex junctions formed by 
trapping a C6o molecule between gold electrodes. In a pre¬ 
vious paper 87 we have demonstrated both experimentally and 
theoretically that Cqq -based nanojunctions show promisingly - 
high values for the thermopower and figure of merit. However 
the temperature dependence of these values and the fluctua¬ 
tions caused by the exact geometrical details have not been 
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FIG. 38: Geometry of a junction having gold (111) electrodes, which 
sandwich a Cqq molecule. The electrodes are terminated as either a 
flat surface (upper panel) or a pyramid (lower panel), and can be 
tilted an angle v relative to each other. 

thoroughly investigated, partly due to the computationally- 
expensive nature of the calculations. Here we show that the 
fast and efficient implementation of GOLLUM allows us to 
undertake a more complete exploration of the transport prop¬ 
erties of these C6o-based junctions. 

The systems of interest consist of two (111) gold leads that 
can be tilted an angle v relative to each other. Each lead is 
terminated using either a flat surface or a pyramid, as shown 
in Fig. ( [38] ). We have performed DFT calculations using the 
code SIESTA, with a double-zeta-polarized basis set, and the 
EDA functional. We have relaxed the molecular geometries 
using a force tolerance of 0.02 eV/A and have found the equi¬ 
librium distance between the leads and the C6o molecule to be 
of about 0.22 nm, depending on the exact orientation of the 
molecule. This result is in good agreement with other previ¬ 
ously reported distance® We have kept this distance fixed 
in all subsequent transport calculations. However, we have 
taken for completeness five possible orientations of the Cqq 
molecule relative to the electrodes. These are: (a) a C-C bond 
between a hexagon and a pentagon facing the Au surface; (b) 
a hexagon facing the Au surface; (c) a pentagon facing the 
Au surface; (d) a bond between two hexagons facing the Au 
surface; and (e) a single atom facing the Au surface. We have 
also tilted one of the electrodes in steps of 15 degrees between 
0 and 60 degrees to see the interference effects caused by the 
exact position of the tip on the surface of the fullerene, recal¬ 
culating the thermoelectric properties at each step. The start¬ 
ing position (for v = 0) for the C6o against the electrodes is 
such that one of its pentagons is facing the Au surfaces. 

The results of our calculations are shown in Figure ( [39] ). By 
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FIG. 39: (Color online) The panels show our thermoelectric re¬ 
sults for the junction in Fig. ( [38] ) having flat or pyramid-terminated 
electrodes; v is the rotation angle. Figs, (a)-(b)-(c)-(d) show two- 
dimensional contour plots of the conductance G (measured in units 
of Go), the thermopower S (measured in /.iV/K ), the thermal con¬ 
ductance k (measured in nW / K), and the figure of merit ZT (di¬ 
mensionless). The vertical and horizontal axes are the temperature 
T measured in Kelvin and the tilting angle v for flat-electrode junc¬ 
tions. Figs, (e)-(f)-(g)-(h) show the same magnitudes for pyramid- 
terminated electrodes. Note that the color code is different for each 
figure. 


taking horizontal cuts through these surfaces, we can see clear 
evidence of quantum interference as the angle changes and the 
tip is repositioned around the fullerene surface. For the con¬ 
ductance G, these oscillations are almost temperature inde¬ 
pendent, whereas in the case of the thermopower S , the ther¬ 
mal conductance k and the electronic figure of merit ZT, these 
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oscillations are almost negligible at low temperatures and then 
grow with T. Comparing the conductance obtained with flat 
electrodes against that with pyramid-terminated electrodes, 
we can clearly see that G decreases substantially when us¬ 
ing the pyramid-like electrodes, in agreement with Ref. ([ 89 ]). 
Furthermore in the case of flat electrodes the conductance 
fluctuates with the angle v by about half an order of magni¬ 
tude, whereas for the pyramid-terminated electrodes we find 
a larger change of almost one order of magnitude. This again 
demonstrates that the pyramid-terminated electrode scans the 
molecular surface like an STM tip, with improved detail, 
while part of these features are blurred when using a flat elec¬ 
trode. The flat-electrode junctions possess conductances val¬ 
ues of about 0.6-1.2 GO, while junctions with pyramid tips 
have conductances of order 0.015-0.15 GO. We overestimate 
the experimental values for G 87 9 321 by about one order of 
magnitude, a known problem associated with the underesti¬ 
mation of the HOMO-LUMO gap inherent to the plain DFT 
approach. As expected, the thermopower S is more sensitive 
to the angle v when using pyramid-terminated junctions com¬ 
pared with the case of flat-surfaced junctions. Interestingly, S 
is quite high especially at the higher temperatures, achieving 
values of of about 100 to 200 /iky iT. 


L. Multi-terminal calculations 

Ab-initio force-relaxation simulations show that it is 
possible to sculpt complex three-dimensional structures of 
nanoscale-scale dimensions by cutting shapes into a graphene 
bilayei®. We find that the edges of the two graphene sheets 
coalesce in order to saturate dangling bonds and to maxi¬ 
mize the degree of sp 2 hybridization. For example, by cut¬ 
ting a cross shape in a bilayer graphene sheet, the resulting 
sculpturene is the three-dimensional crossbar carbon nanotube 
(CNT) shown in Fig. ( [40]) , which is an example of a four- 
terminal electronic device. This four-terminal device is com¬ 
posed of two armchair and two zigzag CNT electrodes. These 
become perfectly periodic CNT leads (shown with blue) far 
enough from the junction. To perform a GOLLUM-based 
four-terminal calculation, we obtain the mean-field Hamilto¬ 
nian of this structure using the SIESTA code, with a double- 
zeta-polarized basis set and a GGA functional 76 . We start with 
the referred cross-shaped bilayer graphene sheet and after re¬ 
laxing the inter-atomic forces to tolerances below 0.02 eV/A, 
we find the crossbar shaped device shown in Figure ( [40] ). 

We feed the resulting ab-initio Hamiltonian into GOLLUM, 
and compute the transmission coefficients T tJ between every 
possible combination of pairs of leads. Notice that the arm¬ 
chair CNT leads are semiconducting, while the zigzag CNTs 
are metallic. We therefore expect different qualitative behav¬ 
iors for the transmission properties among the different arms. 
This is shown in Fig. pT) , where the transmission coefficients 
between the two armchair arms T12 are much smaller than 
those connecting the zigzag arms T 34 . In addition, the figure 
indicates that the central cross area, where the two arms join 
together is not transparent, but introduces strong scattering. 
Similarly, the transmission from lead 1 to lead 3 also shows a 



FIG. 40: (Color online) Four-probe cross bar carbon nanotube de¬ 
vice. Here, LI and L2 are armchair CNTs, while L3 and L4 are 
zigzag CNTs, both with diameters of 5.65 A. The distance between 
the edges of LI and L2 and of L3 and L4 are 50.51 and 49.87 A, 
respectively. 



FIG. 41: (Color online) Transmission coefficients between leads 1 
and 2 T 12 and leads 3 and 4 T 34 (solid black and dashed red lines 
respectively). Energies are measured in eV and referred to the Fermi 
energy of the Scattering Region. 


reduced transmission at low energies due to the semiconduct¬ 
ing behavior of the armchair lead 1, see fig ©• This figure 
shows that T13 = T42 due to the junction symmetry. 


M. Environmental effects on quantum transport 

In the literature, most theoretical analyses of phase- 
coherent transport properties assume that the junction is im¬ 
mersed in vacuum and therefore ignore the effects of the sur¬ 
rounding environment. In contrast, many experiments are car¬ 
ried out under ambient conditions, which can have a marked 
effect on transport properties^ ' If surrounding environmental 
molecules possess a dipole moment, then the scattering re- 
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FIG. 42: (Color online) Transmission coefficients between lead 4 
and lead 2 (T42, solid black line) and between lead 1 and lead 3 (T13, 
dashed red line). Energies are measured in eV and referred to the 
Fermi energy of the Scattering Region. 


gion will be subject to a fluctuating electrostatic field. Pre¬ 
vious work to investigate the impact of environmental water 
on the transport properties of a a single-molecule junction 94 
also took into account the effect of a solvation shell of wa¬ 
ter molecules surrounding the junction. GOLLUM describes 
these effects systematically, by noting that for nanostructures 
such as single-molecule junctions, the timescale for such fluc¬ 
tuations is typically longer than the time taken for an elec¬ 
tron to pass through the device and therefore one can adopt 
the Born-Oppenheimer approximation and freeze the environ¬ 
ment during each electron transit. However, successive elec¬ 
trons experience different environmental snapshots and there¬ 
fore are subjected to different instantaneous mean field Hamil¬ 
tonians leading to different instantaneous conductances. The 
measured (time-averaged) electrical conductance will hence 
be an ensemble average over these snapshots. To obtain a 
series of environmental snapshots, GOLLUM uses classical 
molecular dynamics to describe the environmental molecules 
and for each snapshot, feeds the resulting geometries into 
a DFT code to compute the corresponding self-consistent 
Hamiltonian. The resulting mean field Hamiltonian is then 
used to compute the scattering matrix and related transport 
properties. 

To illustrate this approach, we compute here the ensemble 
averaged conductance of the junction shown in Fig. 
where two pi-stacked monothiol terminated oligophenyle- 
neethynylenes form a bridge between two gold (001) elec¬ 
trodes. The bridging molecule is surrounded by two differ¬ 
ent solvents: decane and 1,4-dioxane, 1,2,4-trichlorobenzene 
(TCB). An example of a junction surrounded by a shell of 
TCB molecules is shown in the lower panel of Fig. ©• 
We have tested here the classical molecular dynamics pack- 
ages LAMMP^l and DLPOL\®, but GOLLUM is flexible 
enough to accept coordinates from other classical Molecu¬ 
lar Dynamics packages. In what follows, we show results 
obtained with LAMMPS, where we have used the Dreiding 
force field to describe the intra- and inter-molecular interac- 




FIG. 43: (Colour online) (top) Geometry of a pi-stacked molecule 
connected to gold electrodes, (bottom) Single snapshot of a MD 
calculation where the molecule is surrounded by TCB solvent 
molecules. 



FIG. 44: (Colour online) Transmission curves of the junction shown 
in Fig. ( [43] ) for 100 snapshots. 


tions and have employed the REAXFF forcefield to obtain the 
initial charges. To create the environment, we place two hun¬ 
dred solvent molecules surrounding the backbone molecule. 
We perform the simulations using a constant temperature and 
volume (NVT) ensemble and subsequently a constant tem¬ 
perature and pressure (NPT) thermostat. We equilibrate the 
junction for 150 ps with 0.1 fs time steps continuously raising 
the temperature to 290 K. We do not include the gold elec¬ 
trodes in the molecular dynamics simulation, so to simulate 
the binding of the anchor groups, we hold the positions of the 
two terminating atoms which connect to the electrodes fixed. 
We record between 350 and 500 snapshots of the junction, 
that have been taken every 2ps. For each snapshot, we feed 
the atomic coordinates into the DFT code SIESTA and gen¬ 
erate the DFT Hamiltonian. We then feed the Hamiltonians 
into the transport code to compute the electrical conductance. 
Some example transmission curves for 100 snapshots of the 
junction with a TCB solvent are shown in the right panel in 
Fig. (43). We note that the the room-temperature dynam¬ 
ics of the atoms at the junction lead to a large spread in the 
transmission curves, and therefore to many different values 
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FIG. 45: (Colour online) Conductance histograms of the junction 
shown in Fig. ( |43] > for two different solvents: decane (left) and TCB 
(right). 



FIG. 46: (Colour online) sketch of a possible CNT based DNA sen¬ 
sor: a piece of a ssDNA is being translocated through a torus-shaped 
sculpturene. 


of the low-voltage conductance. We therefore assemble con¬ 
ductance histograms to help identify the most probable con¬ 
ductance values. The resulting histograms in the presence of 
decane and TCB solvents are shown in Fig. ©• The fact that 
the most probable conductance values are different shows that 
ambient-conditions or liquid-immersed molecular electronics 
experiments are affected by the surrounding solvent. Notice 
that in these simulations we kept fixed the molecule-electrode 
geometry, so the spread in G is due entirely to environmental 
effects. 


N. Nanopore-based DNA nucleobase sensing 





(a) AAT 


(d) GCG 


(b) ATT 



(e) TCG 


(c) CGC 


(f) TTC 




The fact that the transport properties of nanoscale junctions 
depends on the surrounding environment leads to a wide range 
of possible sensing applications. In this section, we demon¬ 
strate the versatility of GOLLUM by showing how it can be 
used to predict the change in conductance of a nanopore, when 
a single DNA strand is trans-located through it. Deoxyri¬ 
bonucleic acid (DNA) is a molecule that encodes the genetic 
instructions used in the development and functioning of all 
known living organisms and many viruses. DNA molecules 
are double-stranded helices, consisting of two long biopoly¬ 
mers composed of simpler units called nucleotides. Each nu¬ 
cleotide is composed of one of the four nucleobases guanine 
(G), adenine (A), thymine (T) and cytosine (C), which are 
attached to a backbone made of alternating sugar and phos¬ 
phate groups. The two polymer strands are bound together 
by non-covalent bonds that link base pairs and are easily sep¬ 
arated to form two single-stranded DNA molecules (ssDNA) 
molecules. DNA sequencing aims at identifying the sequence 
of the DNA bases in a sample of ssDNA. 

Many researchers are actively seeking new methods to se¬ 
quence DNA with improved reliability and scalability and that 
are economically viable. Biological nanopores made from 
protein such as a-hemolysin have been shown experimentally 
to sense the presence of DNAP®, but are also very sensitive 
to temperature and pH, and can only be used within a limited 
voltage bias window". As an alternative, solid state devices 


FIG. 47: (Color online) The six nucleotide sequences of three base 
pairs joined by a DNA backbone, discussed in the text: (a) AAT; (b) 
ATT; (c) CGC; (d) GCG; (e) TCG and (f) TTC. 


which can be integrated into existing semiconducting circuitry 
technology and that are robust to the chemical environment 
have been proposed as sensors" 102 . 

To demonstrate the versatility of GOLLUM, we examine 
here the potential for DNA nucleobase sensing of the sculp¬ 
turene device shown in Lig. (46), which comprises a torus-like 
nanopore connected to two CNT electrode The torus in the 
figure has an inner pore with a diameter of 1.6 nm, whereas 
the leads are two (6,6) armchair nanotubes having a diameter 
of about 5 A. We have selected for our study six short strands 
of ssDNA containing three bases, that are shown in Lig. ©• 
We have first relaxed the coordinates of the nucleotides that 
are threading the pore, using the DLT code SIESTA with a 
double-zeta basis set and a LDA functional. Since the pore 
diameter is slightly larger than the strand width, the strand 
and its nucleotides can adopt different conformations and ori¬ 
entations inside the pore. We accumulate snapshots of these 
different conformations and orientations for each of the six 
ssDNA strands as they trans-locate the pore. Lor each snap¬ 
shot, we compute the current-voltage curve and subtract the 
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V (volts) 


FIG. 48: (Color online) I — V curves for each of the six sequences 
shown in Fig. where the current has been averaged over differ¬ 
ent pore-nucleotide relative angles and the current of the empty pore 
has been subtracted. 


current for the empty pore A/ = I(V) — Iq(V). The cur¬ 
rent averaged over snapshots A (I) for each ssDNA strand is 
plotted in Fig. ( |48| ). The sizable height of the curves demon¬ 
strate that the conductance of the pore is sensitive to the gating 
effect produced by the presence of ssDNA strands inside the 
pore. Furthermore, the different behavior of the curves means 
that, armed with a proper statistical analysis, the sculpturene 
device can distinguish different nucleotide sequences, so that 
this kind of device could be utilized potentially as a discrimi¬ 
nating DNA sensor. 


O. Theoretical simulation of the pulling curves and 
histograms of break-junction experiments. 


A large body of experiments in single-molecule electronics 
is performed using the mechanically-controlled break junction 
(MCBJ) technique, in which a metallic strip is pulled slowly 
until it breaks into two separate pieces. This process enables 
the formation of electrodes with molecular-scale gaps, which 
can be bridged by a single molecule. Experimentally, these 
two electrodes are repeatedly pulled away or pushed towards 
each other. By applying a small bias voltage and recording the 
current passing through the junction, the low-voltage conduc¬ 
tance can be measured as a function of the distance between 
the electrodes. When the distance is small enough, a single 
molecule can bridge the gap between the electrodes and its 
conductance can be measured. In the literature, most theo¬ 
retical studies are confined to small numbers of ideal geome¬ 
tries and binding configurations. In this section, our aim is 
to demonstrate that the versatility of GOLLUM allows us to 
compute whole ’pulling curves’ of conductance versus elec¬ 
trode separation. 

The single-molecule junction that we discuss here is shown 
in Fig. ( |49| ) and consists of gold (111) electrodes. The 



FIG. 49: (Colour online) A snapshot of an opening cycle in a gold- 
bipyridine MCBJ simulation. 



I- H- III. iv. 

FIG. 50: (Color online) Relaxed configurations of the junction shown 
in Fig. ( [49] ) for four different distances d. 


electrodes in the simulation are terminated by pyramids and 
bridged by a bipyridine molecule. To simulate a stretching 
process we have created one hundred geometries of the junc¬ 
tion, each with a different distance d between the center of 
the end atoms of the two leads. To optimize the atomic ar¬ 
rangement of each of the hundred geometries, we start from 
an idealized setup consisting of the two pyramids surrounded 
by vacuum (e.g.: not attached to the gold leads). We place the 
molecule slightly shifted to one side to break the symmetry 
and keep an Au-N bond-length of about 2 A. We then relax 
the inter-atomic forces with the SIESTA code using a GGA 
functional 76 and a double-zeta-polarized basis set until each 
individual force is smaller than 0.02 eV/A. We keep fixed the 
atomic positions of the bottom two layers of the pyramids dur¬ 
ing this geometry optimization. Fig. ( [50] ) shows four of the 
hundred relaxed configurations achieved. 

We then reattach the crystalline gold leads, and impose 
periodic boundary conditions along the plane perpendicu¬ 
lar to the transport direction. We use a SZ basis for the 
gold atoms of the leads, together with a simplified pseudo¬ 
potential, where only the 6s channel is included to speed up 
the simulations. However, we use a double-zeta-polarized ba¬ 
sis set for the atoms at the gold pyramids and in the molecule. 
SIESTA then creates the Hamiltonian of each of the hundred 
junctions that we feed into GOLLUM. 

Figure © shows the low-voltage conductance G ver¬ 
sus the electrode separation d. This ’pulling curve’ shows 
that during the pulling process the conductance possesses a 
plateau, in agreement with many experiments using MCBJs. 
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FIG. 52: (Color online) The device geometry of a double-wall nano¬ 
electromechanical quantum pump. An outer carbon nanotube of 
length L « 50 A surrounds concentrically an inner tube, with an 
inter-layer spacing W ^ 3.4 A corresponding to the van der Waals 
distance. The inner wall remains fixed, while the outer tube is rotated 
about the tube axis by the angle p. A slow variation of p results in 
a parametric current d^Q. Finite charge can be pumped in one rota¬ 
tional cycle depending on the chiralities of the constituent tubes. 


FIG. 51: (Color online) Conductance of the junction shown in Fig. 
( |49| . G is measured in units of Go and plotted as a function of the 
distance d. Each of the conductance points displayed correspond to 
one of the hundred relaxed configurations of the MCBJ simulation. 
The figure also shows the four different distances d corresponding to 
the relaxed geometries displayed in Fig. (50|. 


This simulation also reveals that the aromatic rings contact di¬ 
rectly the gold surface, therefore increasing the molecule-gold 
coupling and the molecular conductance. 


P. Quantum Pumping in Carbon Nanotube Archimedes 
Screws 

So far, all calculated quantities have been obtained from 
the modulus squared of the scattering matrix elements. To 
demonstrate that GOLLUM also provides information about 
transport properties associated with the phases of the scat¬ 
tering matrix, we now examine an example of a quantum 
pump. Quantum pumps are time-dependent electron scatter¬ 
ed, which are able to transport electrons between two exter¬ 
nal reservoirs subjected to the same chemical potential. The 
pump process is adiabatic if the frequency of the pump cycle is 
smaller than the inverse of the characteristic timescale of the 
scatterer, the Wigner delay timd®. Experimenta l 104 * 105 1 and 
theoretical 106 108 studies of adiabatic quantum pumps have ex¬ 
amined the conditions for optimal pumping and the effects of 
noise and dissipation. 

Adiabatic pumping can be understood in terms of the para¬ 
metric deriva tive of the full scattering matrix S at fixed chem¬ 
ical potential 109110 . An adiabatically-slowly time-varying 
scatterer connected by ideal channels to external reservoirs, 
produces a current 

d t Q j (t) = ^£ jj (t,E F ) (95) 

pumped into the j th channel, where £jj is the energy shift 
matrix defined by 

£ (t, e f ) = ih — s* C t, E f ) , (96) 


with S being the full scattering matrix and E F the Fermi en¬ 
ergy. Pumping can occur if the 5-matrix depends on time 
through a parameter p(t). Hence currents can be expressed 
in terms of parametric derivatives 


d t Qj = d<pQj(ip) d t <p(t) (97) 


where d^Qj is the parametric current entering channel j. 

Since GOLLUM gives us access to the full scattering ma¬ 
trix <S, it offers the possibility of investigating adiabatic pump¬ 
ing in nanostructures. We demonstrate this capability by cal¬ 
culating the charge pumped in a double-walled carbon nan¬ 
otube nano-electromechanical device shown in Fig. @nu, 
that mimics the experimental setup of Ref. (1112b . Since an 
electron current travelling along the inner tube can cause a 
chiral outer tube to rotate 113 , the quantum pump shown in 
Fig. (52) represents the inverse effect, in which rotation of 
the outer tube causes a current to flow along the inner tube. 
The position and orientation of the inner tube is kept fixed, 
while the shorter outer tube rotates slowly. The angle p de¬ 
scribes the real space rotation angle of the outer tube and also 
plays the role of the pumping parameter in this system. 

To reveal the rich behavior of this family of quantum 
pumps, Fig. ( [53] ) shows the parametric current, as a function 
of the rotational angle (p for a typical device. Depending on 
the particular angle, charge may be pumped either from left 
to right or vice versa. The integral of this parametric emis- 
sivity within a full parametric cycle of 360° is the number of 
electrons pumped per cycle. 

In Fig. ( [54] ) we show the charge pumped in a (5,5) carbon 
nanotube with a (14,6) outer nanotube rotating slowly about it. 
The average pumped charge clearly drops by several orders of 
magnitude as the Fermi energy is increased from zero. There¬ 
fore for a most efficient pumping, the Fermi energy should 
be close to the Dirac point. Note however, that the pumped 
charge could again increase if the Fermi level is large enough 
to open another channel. Beyond this average behavior, there 
exist numerous sharp peaks in the pumped charge. The lo¬ 
cation of these peaks correlates with Fabry-Perot resonances 
in the reflection coefficient. This suggests that the largest 
pumping occurs at those resonances. In other words, when 
the transmission is high, pumping is low, and vice versa. 
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FIG. 55: A two-terminal tight-binding system defined on a square 
lattice, comprising two leads connected to a disordered scattering 
region. 


FIG. 53: (Color online) Contour plots providing the parametric cur¬ 
rent d^Q for a device consisting of a (5,5) carbon nanotube sur¬ 
rounded by a (14,6) outer nanotube. The outer tube rotates slowly 
around the inner one. The contour plot shows the current as a func¬ 
tion of the rotational angle (p measured in degrees and the Fermi en¬ 
ergy Ep measured in units of the hopping integral 7 between carbon 
atoms. The orange color indicates charge being pumped from left to 
right and vice versa. 



FIG. 54: (Color online) The calculated charge pumped per cycle 
Q P (e) through the device discussed in Fig. ( [53] ) as a function of the 
Fermi energy measured in units of the hopping integral 7 between 
carbon atoms at the CNTs. The green solid line shows the charge 
pumped towards the left; the red solid line shown the charge pumped 
towards the right. At certain energies, the pumped charge is very 
high. These peaks correlate with the Fabry-Perot resonances in the 
reflection coefficient R. R is shown as a solid blue line. 



FIG. 56: Conductivity a and transmission coefficient T(Ef) as a 
function of the number L of unit cells in the transport direction, for 
the two probe square lattice shown in the previous figure. 


nian of the system has a single orbital per site, with nearest 
neighbour couplings. 7 = — 1 eV. The site energies within the 
leads are £$ = 0 eV, while the random site energies within the 
scattering region are uniformly distributed over the interval 


[-1.6,1.6] eV. 

Figure ( [56] ) shows the the ensemble-averaged conductiv¬ 
ity (a) and transmission coefficient ( T(Ep )) for the system 
shown in Fig. (55). It contains three regions. Within the bal¬ 
listic regime between L = 0 and approximately L = 20, the 
conductivity increases linearly with length. In the diffusive re¬ 
gion (L = 40 — 80), the conductivity exhibits ohmic behavior 
and is almost independent of length. Finally for L greater than 
100, there is a cross over to the Anderson localized regime. 


Q. Transport in disordered systems: ballistic, diffusive and 
localized behavior 

Finally, to demonstrate that GOLLUM can handle the dis¬ 
ordered systems, we calculate the ensemble-averaged conduc¬ 
tivity a of a two-terminal system on a square lattice, with 
leads attached to a disordered scattering region as shown in 
Fig. ( [55] ). The width of the system is W = 11 unit cells and 
the length is varied between L = 1 to L = 500 unit cells. The 
conductivity is defined as a = T(Ep)W/L , where T(Ep) is 
the transmission from the left lead to right lead evaluated at 
the Fermi energy, Ep = 0.5 eV. The tight-binding Hamilto- 


R. Impact of the spin-orbit interaction in the transport 
properties of nickel chains 

We end this article by showing how the spin-orbit inter¬ 
action induces gaps at certain band crossings in the one¬ 
dimensional electronic structure of infinite nickel chains. 
These gaps may appear or not depending on the orientation 
of the atomic spins relative to the axis of the chain. They lead 
to dips in the transmission coefficient T(E) of the chain at the 
gap energies. 

We have simulated linear nickel chains using the DFT pro- 
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will soon be freely available from the following web site 
http://www.physics.lancs.ac.uk/gollum and the authors of this 
article are available to help potential users access the code. 
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Appendix A: Procedures used to regularize K i 


FIG. 57: Transmission coefficient of a linear nickel chain oriented 
along the z-axis. The solid black line shows T(E) when the atomic 
spins are oriented along the chains axis. The dashed red line shows 
T(E) when the atomic spins are oriented in the plane perpendicular 
to the chain axis. A blue dotted line showing T(E) when the spin- 
orbit interaction is switched off falls on top of the black line. 


gram SIESTA. The chains have a single atom per unit cell and 
are oriented along the z-axis. We have used a standard set of 
pseudopotential parameters, a double-zeta basis set and sim¬ 
ple LDA for the exchange-correlation potential. 

We have checked that the electronic structure and T(E) are 
the same for any spin orientation if the spin-orbit interaction 
is set to zero, as is should due rotational invariance. How¬ 
ever, if the spin-orbit interaction is switched on, then a finite 
yet small magnetic anisotropy barrier appears. We have found 
that if we choose the atomic spins to lie along the chain axis, 
then there are no spin-orbit gaps close to the Fermi energy. As 
a consequence, the transmission coefficients with and with¬ 
out the spin-orbit interaction are indistinguishable from each 
other (as shown in Fig. ([57])). In contrast, when the atomic 
spins are oriented in a plane perpendicular to the chain axis, 
then several small gaps open around the Fermi energy. These 
gaps are seen as dips in Fig. ([57 ). 


IV. CONCLUSION 

We have developed a new quantum transport code, which is 
fast, easy to use and versatile. This flexibility has been demon¬ 
strated by presenting a wide range of example calculations, 
encompassing charge, spin and thermal transport, corrections 
to density functional theory such as LDA+U and spectral ad¬ 
justments, transport in the presence of non-collinear mag¬ 
netism, the quantum-Hall effect, Kondo and Coulomb block¬ 
ade effects, finite-voltage transport, multi-terminal trans¬ 
port, quantum pumps, superconducting nanostructures, en¬ 
vironmental effects and pulling curves and conductance his¬ 
tograms for mechanically-controlled-break-junction experi¬ 
ments. Further developments are in the pipeline, includ¬ 
ing the incorporation of phonon transport. GOLLUM 


We have described in section II.A.2 the method employed 
by GOLLUM to find the surface Green’s function Gs of each 
lead. However, the solution of Eq. © gives with some 
frequency numerical inaccuracies which render the method 
useless as it stands. These inaccuracies are caused by the 
highly non-singular behavior of the Hamiltonian matrix K\ 
connecting adjacent PLs. We discuss here the adaptation of 
the method described in Refs. (19137b that GOLLUM uses to 
regularize K\. Mathematically, we perform an SVD decom¬ 
position of this N x N matrix, 

K 1 = USV f (Al) 

where U and V are unitary matrix and S' is a diagonal matrix 
containing the eigenvalues A of Numerical algorithms 
usually arrange them in descending order. The condition num¬ 
ber of K\ is defined as the ratio k = A max /A min between 
the maximum and minimum eigenvalues of K\. k determines 
how singular is K f 1 and therefore the propensity to suffer 
inaccuracies when handling K\. Small eigenvalues A appear 
whenever K\ is very sparse. Physically, this is originated for 
example if the PL are very long so that a large fraction of hop¬ 
ping integrals (matrix elements of K\) is zero. Our first pro¬ 
cedure to regularize K\ consists in adding a real or complex 
random matrix to K\. We have found that this is frequently 
enough to render a regular K\ matrix. If this first procedure 
fails this is because the orbitals involved do not play a role in 
the transport properties of the lead and should be decimated 
out, so that the dimensions of the K\ matrix are reduced. 

Notice that reducing the dimensions of K\ has the advan¬ 
tage that the computation of the surface Green’s function Gs 
is much lighter. However, the procedure must be performed 
with some care as finally Gs must connect with the corre¬ 
sponding TPL of the EM branch, whose matrices have dimen¬ 
sions N x N. Explicitly, we write the Schroedinger equation 
of the infinite chain of the corresponding Lead: 

K 0 C n + K\ C n + 1 + K i C n - 1 = 0 (A2) 

where C n = e lkna C(k ), n labels the PL and runs from — oo 
to +oo and we assume that the chain will be chopped off at 
the n o PL and then connected to the TPL of the EM. We first 
perform the SVD decomposition of K\ described in Eq. 
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We then set to zero all eigenvalues smaller than a given toler- We transform Eq. ( |A2| ) for all sites up to site no — 1 as follows: 

ance tol. We have checked that setting tol = 10“ 8 — 10 -9 

provides unproblematic K\ and K_ i = K\. Let us assume 

that we set to zero D eigenvalues of Ki, so that M = N — D 

remain non-zero. We now construct the matrices 


K, = U I K ' M ° I V' 

0 0 


I<[ = (Kirf = F f KiV = F f F S = 


Ki ,m 0 
P 0 


K 


A'i = 


CL = 


F f K 0 V = 
k\ = F f i<i = 


Ko M W 


Wt Ko,d 
Vi f 2 


F 3 V 4 


A 

B 


V* Cn = | C ' M ' n 

C'd.u 


(A3) 


Ft (K 0 C n + K x C n+1 + K-i) V Ft c n _! = Ki C' n + I<[ C’ n+l + K’_ x C' n _ x = 0 (A4) 


while the equation for sites no and no + 1 will be We now decimate out C' D up to site no, arriving to the new set 

of equations 

K c' n o + kl i c; 0 _i + k x c n0+ 1 = o 

K 0 C n o + i + K- 1 C no + Ki C n +2 = 0 (A5) 


{Kom ~ (F f K-d IF - Ft K ~1 d P ) C ' Mtn + (K 1M ~ IFt p) C^ n+1 + (< M - pt A' ( ^ If) = 0 n 

-V- V* V* 

K 0 new iV“ ew 

(A>,m - IFt k -1 d w _ pi k -1 d P ) C M nQ + (a - IFt p) c; 0+1 + (kI m - pt K-i, if) = 0 

(p 0 - pt k 0 , d b ) c; 0+1 + (At - pt Kq d W ) + Kl c no+2 = o 

This set of equations generates the new Hamiltonians of each lead, and connects it with the extended molecule. Overlap matrices 
must also be decimated as they enter into the expression for the group velocity, Eq. ©• The denominator in this equation looks 
like 


C(k )t (S 0 + Si e ika + S_i e ~ ika ) C(k) = C'(fc)' t ( S' 0 + S[ e ika + S'_! e~ ika ) C(k)' 

where we have applied the SVD transformation to the overlap matrices 

Si = F t S' 0 F = 


Si = F t S'iF = 


(A7) 



D 0 
F 0 


(A8) 


< n 0 


(A6) 
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By decimating out the unwanted degrees of freedom Cd, we arrive at the following expressions 

SS ew = A-W' K~b B + B t K-}, W + K~b \,C K~bW + pt K~b C K~b P ~ P ] K~b E - E* A' ( ^ 
5i new = D - B' K-b P - W/t K~b E + W/t K-^CK~b P 


P 

(A9) 
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